Ultrahigh Energy Cosmic Rays: Facts, Myths, and Legends

Ultrahigh Energy Cosmic Rays: Facts, Myths, and Legends

Luis Alfredo Anchordoqui

University of Wisconsin-Milwaukee, USA

Abstract

This is a written version of a series of lectures aimed at graduate students in
astrophysics/particle theory/particle experiment. In the first part, we explain
the important progress made in recent years towards understanding the ex-
perimental data on cosmic rays with energies > 10 8 GeV. We begin with a
brief survey of the available data, including a description of the energy spec-
trum, mass composition, and arrival directions. At this point we also give a
short overview of experimental techniques. After that, we introduce the fun-
damentals of acceleration and propagation in order to discuss the conjectured
nearby cosmic ray sources, and emphasize some of the prospects for a new
(multi-particle) astronomy. Next, we survey the state of the art regarding the
ultrahigh energy cosmic neutrinos which should be produced in association
with the observed cosmic rays. In the second part, we summarize the phe-
nomenology of cosmic ray air showers. We explain the hadronic interaction
models used to extrapolate results from collider data to ultrahigh energies, and
describe the prospects for insights into forward physics at the Large Hadron
Collider (LHC). We also explain the main electromagnetic processes that gov-
ern the longitudinal shower evolution. Armed with these two principal shower
ingredients and motivation from the underlying physics, we describe the dif-
ferent methods proposed to distinguish primary species. In the last part, we
outline how ultrahigh energy cosmic ray interactions can be used to probe new
physics beyond the electroweak scaled

^Lecture notes from the 6th CERN-Latin- American School of High-Energy Physics, Natal, Brazil, March - April, 2011.
http : //physicschool . web . cern . ch/PhysicSchool/CLASHEP/CLASHEP201 1/.

Contents

1 Setting the stage 2

2 Frontiers of multi-messenger astronomy 3

2.1 A century of cosmic ray observations 3

2.2 Origin of ultrahigh energy cosmic rays 16

2.3 Energy losses of baryonic cosmic rays on the pervasive radiation fields 24

2.4 Diffuse propagation of protons in a magnetized Local Supercluster 30

2.5 Ultrahigh energy cosmic neutrinos 35

3 Phenomenology of extensive air showers 50

3.1 Systematic uncertainties in air shower measurements from hadronic interaction models 50

3.2 Electromagnetic processes 61

3.3 Paper- and-pencil air shower modeling 65

4 Searches for new physics beyond the electroweak scale at ~ 250 TeV 69

4.1 General idea 69

4.2 v acceptance and systematic uncertainties 71

4.3 Auger discovery reach 72

A Cosmogenic /3-DK and A* processes 76

B Energy density of the EM cascade 78

C TOTal Elastic and diffractive cross section Measurement 80

D Raiders of the lost holy grail 81

1 Setting the stage

For biological reasons our perception of the Universe is based on the observation of photons, most

trivially by staring at the night-sky with our bare eyes. Conventional astronomy covers many orders of
magnitude in photon wavelengths, from 10 4 cm radio-waves to 1CT 14 cm gamma rays of GeV energy.

This 60 octave span in photon frequency allows for a dramatic expansion of our observational capacity
beyond the approximately one octave perceivable by the human eye. Photons are not, however, the
only messengers of astrophysical processes; we can also observe cosmic rays and neutrinos (and maybe
also gravitons in the not-so-distant future). Particle astronomy may be feasible for neutral particles or
possibly charged particles with energies high enough to render their trajectories magnetically rigid. This
new astronomy can probe the extreme high energy behavior of distant sources and perhaps provide a
window into optically opaque regions of the Universe. In these lectures we will focus attention on the
highest energy particles ever observed. These ultrahigh energy cosmic rays (UHECRs) carry about seven
orders of magnitude more energy than the LHC beam. We will first discuss the requirements for cosmic
ray acceleration and propagation in the intergalactic space. After that we will discuss the properties of
the particle cascades initiated when UHECRs interact in the atmosphere. Finally we outline strategies
to search for physics beyond the highly successful but conceptually incomplete Standard Model (SM) of
weak, electromagnetic, and strong interactions.

Before proceeding, we pause to present our notation. Unless otherwise stated, we work with
natural (particle physicist's) Heaviside-Lorentz (HL) units with

h = c = k = £Q = HQ = l.

The fine structure constant is a = e 2 / (ineohc) ~ 1/137. All SI units can then be expressed in electron
Volt (eV), namely

1 m ~ 5.1 x 10 6 eV" 1 , 1 s ~ 1.5 x 10 15 eV" 1 , 1 kg ~ 5.6 x 10 35 eV ,

1 A ~ 1244 eV , 1 G ~ 1.95 x 10~ 2 eV 2 , 1 K ~ 8.62 x 10" 5 eV .

2

Table 1: Symbols and units of most common quantities.

symbol quantity unit

J diffuse flux GeV _i cm -2 s _1 sr~

<3> diffuse neutrino flux (upper bound) GeV -1 cm -2 s _1 sr~

(j) point-source neutrino flux (upper bound) GeV -1 cm -2 s _1

L source luminosity ergs -1

C integrated luminosity fb -1 = 10 -39 cm 2

Q source emissivity GeV -1 s _1

Q source emissivity per volume GeV -1 cm -3 s -1

Electromagnetic processes in astrophysical environments are often described in terms of Gauss (G)
units. For a comparison of formulas used in the literature we note the conversion, (e 2 )nL = 47r(e 2 )c,
(£> 2 )hl = (£> 2 )g/47t and (£? 2 )hl = {E 2 )q/Ati. To avoid confusion we will present most of the formu-
las in terms of 'invariant' quantities, i.e. eB, eE and the fine- structure constant a. Distances are generally
measured in Mpc ~ 3.08 x 10 22 m. Extreme energies are sometimes expressed in EeV = 10 18 eV. The
symbols and units of most common quantities are summarized in Table 1.

2 Frontiers of multi-messenger astronomy
2.1 A century of cosmic ray observations

In 1912 Victor Hess carried out a series of pioneering balloon flights during which he measured the
levels of ionizing radiation as high as 5 km above the Earth's surface [1]. His discovery of increased
radiation at high altitude revealed that we are bombarded by ionizing particles from above. These CR
particles are now known to consist primarily of protons, helium, carbon, nitrogen and other heavy ions
up to iron. Measurements of energy and isotropy showed conclusively that one obvious source, the Sun,
is not the main source. Only below 100 MeV kinetic energy or so, where the solar wind shields protons
coming from outside the solar system, does the Sun dominate the observed proton flux. Spacecraft
missions far out into the solar system, well away from the confusing effects of the Earth's atmosphere
and magnetosphere, confirm that the abundances around 1 GeV are strikingly similar to those found in
the ordinary material of the solar system. Exceptions are the overabundance of elements like lithium,
beryllium, and boron, originating from the spallation of heavier nuclei in the interstellar medium.

Above ~ 10 5 GeV, the rate of CR primaries is less than one particle per square meter per year and
direct observation in the upper layers of the atmosphere (balloon or aircraft), or even above (spacecraft)
is inefficient. Only ground-based experiments with large apertures and long exposure times can hope to
acquire a significant number of events. Such experiments exploit the atmosphere as a giant calorimeter.
The incident cosmic radiation interacts with the atomic nuclei of air molecules and produces air showers
which spread out over large areas. Already in 1938, Pierre Auger concluded from the size of the air
showers that the spectrum extends up to and perhaps beyond 10 6 GeV [2, 3]. In recent years, substantial
progress has been made in measuring the extraordinarily low flux at the high end of the spectrum.

There are two primary detection methods that have been successfully used in high exposure exper-
iments. In the following paragraph we provide a terse overview of these approaches. For an authoritative
review on experimental techniques and historical perspective see [4]. For a more recent comprehensive
update including future prospects see [5].

The size of an extensive air shower (EAS) at sea-level depends on the primary energy and arrival
direction. For showers of UHECRs, the cascade is typically several hundreds of meters in diameter
and contains millions of secondary particles. Secondary electrons and muons produced in the decay of
pions may be observed in scintillation counters or alternatively by the Cherenkov light emitted in water
tanks. The separation of these detectors may range from 10 m to 1 km depending on the CR energy

3

Number

TeV muons

Fig. 1: Particles interacting near the top of the atmosphere initiate an electromagnetic and hadronic cascade. Its
profile is shown on the right. The different detection methods are illustrated. Mirrors collect the Cerenkov and
nitrogen fluorescent light, arrays of detectors sample the shower reaching the ground, and underground detectors
identify the muon component of the shower. From Ref . [6] .

and the optimal cost-efficiency of the detection array. The shower core and hence arrival direction can
be estimated by the relative arrival time and density of particles in the grid of detectors. The lateral
particle density of the showers can be used to calibrate the primary energy. Another well-established
method of detection involves measurement of the shower longitudinal development (number of particles
versus atmospheric depth, shown schematically in Fig. 1) by sensing the fluorescence light produced via
interactions of the charged particles in the atmosphere. The emitted light is typically in the 300 - 400 nm
ultraviolet range to which the atmosphere is quite transparent. Under favorable atmospheric conditions,
EASs can be detected at distances as large as 20 km, about 2 attenuation lengths in a standard desert
atmosphere at ground level. However, observations can only be done on clear Moonless nights, resulting
in an average 10% duty cycle.

In these lectures we concentrate on results from the High Resolution Fly's Eye experiment and the
Pierre Auger Observatory to which we will often refer as HiRes and Auger. HiRes [7] was an up-scaled
version of the pioneer Fly's Eye experiment [8]. The facility was comprised of two air fluorescence
detector sites separated by 12.6 km. It was located at the US Army Dugway Proving Ground in the
state of Utah at 40.00° N, 113° W, and atmospheric depth of 870 g/cm 2 . Even though the two detectors
(HiRes-I and HiRes-II) could trigger and reconstruct events independently, the observatory was designed
to measure the fluorescence light stereoscopically. The stereo mode allows accurate reconstruction of the
shower geometry (with precision of 0.4°) and provides valuable information and cross checks about the
atmospheric conditions at the time the showers were detected. HiRes-I and HiRes-II collected data until

4

April 2006 for an accumulated exposure in stereoscopic mode of 3, 460 hours . The monocular mode
had better statistical power and covered a much wider energy range.

The Pierre Auger Observatory [9] is designed to measure the properties of EASs produced by CRs
at the highest energies, above about 10 9 GeV. It features a large aperture to gather a significant sample
of these rare events, as well as complementary detection techniques to mitigate some of the systematic
uncertainties associated with deducing properties of CRs from air shower observables.

The observatory is located on the vast plain of Pampa Amarilla, near the town of Malargue in
Mendoza Province, Argentina (35.1° - 35.5° S, 69.6° W, and atmospheric depth of 875 g/cm 2 ). The
experiment began collecting data in 2004, with construction of the baseline design completed by 2008.
As of October 2010, Auger had collected in excess of 20, 000 km 2 sr yr in exposure, significantly more
exposure than other cosmic ray observatories combined. Two types of instruments are employed. Particle
detectors on the ground sample air shower fronts as they arrive at the Earth's surface, while fluorescence
telescopes measure the light produced as air shower particles excite atmospheric nitrogen.

The surface array [10] comprises 1, 600 surface detector (SD) stations, each consisting of a tank
filled with 12 tons of water and instrumented with 3 photomultiplier tubes which detect the Cherenkov
light produced as particles traverse the water. The signals from the photomultipliers are read out with
flash analog to digital converters at 40 MHz and timestamped by a GPS unit, allowing for detailed study
of the arrival time profile of shower particles. The tanks are arranged on a triangular grid with a 1.5 km
spacing, covering about 3, 000 km 2 . The surface array operates with close to a 100% duty cycle, and the
acceptance for events above 10 9 5 GeV is nearly 100%.

The fluorescence detector (FD) system [11] consists of 4 buildings, each housing 6 telescopes
which overlook the surface array. Each telescope employs an 11 m 2 segmented mirror to focus the
fluorescence light entering through a 2.2 m diaphragm onto a camera which pixelizes the image using
440 photomultiplier tubes. The photomultiplier signals are digitized at 10 MHz, providing a time profile
of the shower as it develops in the atmosphere. The FD can be operated only when the sky is dark and
clear, and has a duty cycle of 10-15%. In contrast to the SD acceptance, the acceptance of FD events
depends strongly on energy [12], extending down to about 10 9 GeV.

The two detector systems provide complementary information, as the SD measures the lateral
distribution and time structure of shower particles arriving at the ground, and the FD measures the lon-
gitudinal development of the shower in the atmosphere. A subset of showers is observed simultaneously
by the SD and FD. These "hybrid" events are very precisely measured and provide an invaluable cali-
bration tool. In particular, the FD allows for a roughly colorimetric measurement of the shower energy
since the amount of fluorescence light generated is proportional to the energy deposited along the shower
path; in contrast, extracting the shower energy via analysis of particle densities at the ground relies on
predictions from hadronic interaction models describing physics at energies beyond those accessible to
current experiments. Hybrid events can therefore be exploited to set a model-independent energy scale
for the SD array, which in turn has access to a greater data sample than the FD due to the greater live
time.

The Extreme Universe Space Observatory (EUSO) is currently being considered by the Japan
Aerospace Exploration Agency for a possible payload on the Japanese Experimental Module (JEM) of
the International Space Station [13]. The mission is currently scheduled for 2017 [14]. The launch
will be provided by the H-II Transfer Vehicle Kounotori. Looking straight down, the UV telescope of
JEM-EUSO will have 1 km 2 resolution on the Earth's surface, which provides an angular resolution of
2.5°. The surface area expected to be covered on Earth is about 160, 000 km 2 , with a duty cicle of order
of 10%. The detector will also operate in a tilted mode that will increase the viewing area by a factor
of up to 5, but decreasing its resolution. The telescope will be equiped with devices that measure the
transparency of the atmosphere and the existence of clouds.

In the remainder of the section, we describe recent results from HiRes and Auger, including the

5

10~

CO

CI

I

10 4 -

>

CD

o

r-

10-*

10

13

Grigorov a
JACEE v
MGU v
TienShan o
Tibet07 •
Akeno o
CASA/MIA eb
Hegra o
Flys Eye •
Agasa ♦
HiResl #
HiRes2 □
Auger SD o
Auger hybrid a
Kascade ffi

10

14

10

15

10 16 10 17
E [eV]

10 J

10 1

10

20

Fig. 2: All-particle spectrum of cosmic rays. From Ref. [15].

measurement of the cosmic ray energy spectrum, composition, and searches for anisotropy in the cosmic
ray arrival directions.

2.1.1 Energy spectrum

The CR spectrum spans over roughly 1 1 decades of energy. Continuously running monitoring using so-
phisticated equipment on high altitude balloons and ingenious installations on the Earth's surface encom-
pass a plummeting flux that goes down from 10 4 m -2 s _1 at ^ 1 GeV to 10 -2 km -2 yr _1 at ~ 10 11 GeV.
Its shape is remarkably featureless, with little deviation from a constant power law (J oc i?~ 7 , with
7^3) across this large energy range. The small changes in the power index, 7, are conveniently visu-
alized taking the product of the flux with some power of the energy. In this case the spectrum reveals a
leg-like structure as it is sketched in Fig. 2. The anatomy of this cosmic leg - its changes in slope, mass
composition or arrival direction - reflects the various aspects of CR propagation, production and source
distribution.

A steepening of the spectrum (7 ~ 2.7 — » 3.1) at an energy of about 10 6 5 GeV is known as the
cosmic ray knee. Composition measurements in cosmic ray observatories indicate that this feature of the
spectrum is composed of the subsequent fall-off of Galactic nuclear components with maximal energy
E/Z [16-18]. This scaling with atomic number Z is expected for particle acceleration in a magnetically
confining environment, which is only effective as long as the particle's Larmor radius is smaller than the
size of the accelerator. If this interpretation holds, the Galactic contribution in CRs can not extend much
further than 10 8 GeV, assuming iron {Z = 26) as the heaviest component. However, the observational
data at these energies is inconclusive and the end-point of Galactic CRs has not been pinned down so far.
For a survey of spectral features at lower energies see [19].

6

The onset of an extragalactic contribution could be signaled by the so-called second knee - a
further steepening (7 ~ 3.1 —> 3.2) of the spectrum at about 10 8 7 GeV. Note that extragalactic CRs are
subject to collisions with the interstellar medium during their propagation over cosmological distances.
Depending on the initial chemical composition, these particle- specific interactions will be imprinted in
the spectrum observed on Earth. It has been argued [20, 21] that an extragalactic proton population
with a simple power-law injection spectrum may reproduce the spectrum above the second knee. In
these models the flattening (7 ~ 3.2 —> 2.7) in the spectrum at around 10 9 5 GeV - the so-called ankle
- can be identified as a "dip" from e + e~ pair production together with a "pile-up" of protons from
pion photoproduction. However, this feature relies on a proton dominance in extragalactic cosmic rays
since heavier nuclei like oxygen or iron have different energy loss properties in the cosmic microwave
background (CMB) and mixed compositions, in general will not reproduce the spectral features [22]. A
cross-over at higher energies is also feasible: above the ankle the Larmor radius of a proton in the galactic
magnetic field exceeds the size of the Galaxy and it is generally assumed that an extragalactic component
dominates the spectrum at these energies [23]. Moreover, the galactic-extragalactic transition ought to be
accompanied by the appearence of spectral features, e.g. two power-law contributions would naturally
produce a flattening in the spectrum if the harder component dominates at lower energies. Hence, the
ankle seems to be the natural candidate for this transition.

Shortly after the CMB was discovered [24], Greisen [25], Zatsepin, and Kuzmin [26] (GZK)
pointed out that the relic photons make the universe opaque to CRs of sufficiently high energy. This
occurs, for example, for protons with energies beyond the photopion production threshold,

where m p (m n ) denotes the proton (pion) mass and c^cmb ~ 10 -3 eV is a typical CMB photon en-
ergy. After pion production, the proton (or perhaps, instead, a neutron) emerges with at least 50% of
the incoming energy. This implies that the nucleon energy changes by an e-folding after a propagation
distance < (a P1 n 7 y^) -1 ~ 15 Mpc. Here, n 7 « 410 cm -3 is the number density of the CMB photons,
a vl > 0.1 mb is the photopion production cross section, and y n is the average energy fraction (in the lab-
oratory system) lost by a nucleon per interaction. For heavy nuclei, the giant dipole resonance (GDR) can
be excited at similar total energies and hence, for example, iron nuclei do not survive fragmentation over
comparable distances. Additionally, the survival probability for extremely high energy 10 11 GeV) 7-
rays (propagating on magnetic fields » 10 -11 G) to a distance d, p(> d) « exp[— d/6.6 Mpc], becomes
less than 10 -4 after traversing a distance of 50 Mpc. All in all, as our horizon shrinks dramatically for
energies > 10 11 GeV, one would expect a sudden suppression in the energy spectrum if the CR sources
follow a cosmological distribution.

At the beginning of summer 2002, in a pioneering paper Bahcall and Waxman [27] noted that
the energy spectra of CRs reported by the AGASA, the Fly's Eye, the Haverah Park, the HiRes, and
the Yakutsk collaborations (see Fig. 3) are consistent with the expected GZK suppression at ~ 3. So-
level according to the Fly's Eye normalization, increasing up to ~ 8a if the selected normalization
is that of Yakutsk. Five years later, the HiRes Collaboration reported a suppression of the CR flux
above E = [5.6 ± 0.5(stat) ± 0.9(syst)] x 10 10 GeV, with 5.3a significance [28]. The spectral index
of the flux steepens from 2.81 ± 0.03 to 5.1 ± 0.7. The discovery of the GZK suppression has been
confirmed by the Pierre Auger Collaboration, measuring 7 = 2.69 ± 0.2(stat) ± 0.06(syst) and 7 =
4.2 ± 0.4(stat) ± 0.06(syst) below and above E = 4.0 x 10 10 GeV, respectively (the systematic
uncertainty in the energy determination is estimated as 22%) [29].

Last year, an updated Auger measurement of the energy spectrum was published [30], corre-
sponding to a surface array exposure of 12, 790 km 2 sr yr. This measurement, combining both hybrid
and SD-only events, is shown in Fig. 4. The ankle feature and flux suppression are clearly visible.
A broken power law fit to the spectrum shows that the break corresponding to the ankle is located at

(1)

10'

26

O Yakutsk T500

□ Yakutsk T1 000

■ Akeno 1 km 2

• AGASA

* I I Haverah Park
x SUGAR
* MSU

>

CD

10'

25

O O

o o

o

CD
CO

w 10'

,24

*****

* HiRes-MIA

a HiRes I mono

▼ HiRes II mono

A Fly's Eye mono

Fly's Eye stereo

10'

23

_i I I i i i i

* Auger 2008

_i i i i i i i i

10 1

VI 9

10'

Energy (eV)

10'

20

Fig. 3: Comparison of flux measurements scaled by E 3 . Only statistical errors are shown. Shown are the data
of AGASA, Akeno, Auger, Fly's Eye, Haverah Park, HiRes-MIA, HiRes Fly's Eye, MSU, SUGAR, and Yakutsk.
Yakutsk T500 (trigger 500) refers to the smaller sub-array of the experiment with 500m detector spacing and
T1000 (trigger 1000) to the array with 1000m detector distance. In case several analyses of the same data set are
available, only the most recent results are included in the plot. The shaded area, depicting the results of the analysis
of the Haverah Park data, accounts for some systematic uncertainties by assuming extreme elemental compositions,
either fully iron or proton dominated. The highest energy point (Fly's Eye monocular observation) corresponds to
the highest energy event. For sake of clarity, upper limits are not shown. The data of the MSU array are included
to show the connection of the high-energy measurements to lower energy data covering the knee of the cosmic-ray
spectrum. FromRef. [19].

log 10 (£7/eV) = 18.61 ± 0.01 with 7 = 3.26 ± 0.04 before the break and 7 = 2.59 ± 0.02 after it. The
break corresponding to the suppression is located at \og 10 (E/eV) = 19.46 ±0.03. Compared to a power
law extrapolation, the significance of the suppression is greater than 20a.

2.1.2 Primary species

Unfortunately, because of the highly indirect method of measurement, extracting precise information
from EASs has proved to be exceedingly difficult. The most fundamental problem is that the first gen-
erations of particles in the cascade are subject to large inherent fluctuations and consequently this limits
the event-by-event energy resolution of the experiments. In addition, the center-of-mass (cm.) energy
of the first few cascade steps is well beyond any reached in collider experiments. Therefore, one needs
to rely on hadronic interaction models that attempt to extrapolate, using different mixtures of theory and
phenomenology, our understanding of particle physics.

The longitudinal development has a well defined maximum, usually referred to as X max , which
increases with primary energy as more cascade generations are required to degrade the secondary particle
energies. Evaluating X max is a fundamental part of many of the composition studies done by detecting
air showers. For showers of a given total energy E, heavier nuclei have smaller X max because the shower
is already subdivided into A nucleons when it enters the atmosphere. The average depth of maximum

8

log io (E/eV)

18 18.5 19 19.5 20 20.5

10 18 10 19 10 20

Energy [eV]

Fig. 4: Combined spectrum from Auger (hybrid and SD events) and the stereo spectrum HiReS. The Auger
systematic uncertainty of the flux scaled by E 3 , due to the uncertainty of the energy scale of 22%, is indicated by
arrows. The results of the two experiments are consistent within systematic uncertainties. From Ref. [30].

(I max ) scales approximately as \n(E/A) [31]. Therefore, since (X max ) can be determined directly from
the longitudinal shower profiles measured with a fluorescence detector, the composition can be extracted
after estimating E from the total fluorescence yield. Indeed, the parameter often measured is Diq, the
rate of change of (X max ) per decade of energy.

Photons penetrate quite deeply into the atmosphere due to decreased secondary multiplicities and
suppression of cross sections by the Landau-Pomeranchuk-Migdal (LPM) effect [32,33]. Indeed, it is
rather easier to distinguish photons from protons and iron than protons and iron are to distinguish from
one another. For example, at 10 10 GeV, the (X max ) for a photon is about 1000 g/cm 2 , while for protons
and iron the numbers are 800 g/cm 2 and 700 g/cm 2 , respectively.

Searches for photon primaries have been conducted using both the surface and fluorescence istru-
ments of Auger. While analysis of the fluorescence data exploits the direct view of shower development,
analysis of data from the surface detector relies on measurement of quantities which are indirectly related
to the X max , such as the signal risetime at 1000 m from the shower core and the curvature of the shower
front. Presently, the 95% CL upper limits on the fraction of CR photons above 2, 3, 5, 10, 20, and
40 x 10 9 GeV are 3.8%, 2.4%, 3.5%, 2.0%, 5.1%, and 31%, respectively. Further details on the analysis
procedures can be found in [34-36]. In Fig. 5 these upper limits are compared with predictions of the
cosmogenic photon flux.

In Fig. 6 we show the variation of (X max ) with primary energy as measured by several experi-
ments. Interpreting the results of these measurements relies on comparison to the predictions of high

9

Fig. 5: Upper limits on the photon fraction reported by the Pierre Auger Collaboration. A prediction for the
cosmogenic photon flux is also shown for comparison (details on the calculation are given Sec. 2.5.2) The color
and line coding is the same as the one in Fig. B.l. This figure is courtesy of Markus Ahlers.

energy hadronic interaction models. As one can see in Fig. 6, there is considerable variation in pre-
dictions among the different interaction models. For 1.6 x 10 9 GeV < E < 6.3 x 10 10 GeV, the
HiRes data are consistent with a constant elongation rate d(X max ) /d(log(E)) = 47.9 ± 6.0(stat) ±
3.2(sys) g/cm 2 /decade [37]. The inference from the HiRes data is therefore a change in cosmic ray
composition, from heavy nuclei to protons, at E ~ 10 8 6 GeV [38]. This is an order of magnitude lower
in energy than the previous crossover deduced from the Fly's Eye data [39]. On the other hand, Auger
measurements, interpreted with current hadronic interaction models, seem to favor a mixed (protons +
nuclei) composition at energies above 10 8 6 GeV [40]. However, uncertainties in the extrapolation of the
proton-air interaction - cross section [41] and elasticity and multiplicity of secondaries [42] - from ac-
celerator measurements to the high energies characteristic for air showers are large enough to undermine
any definite conclusion on the chemical composition.

2.1.3 Distribution of arrival directions

There exists "lore" that convinces us that the highest energy CRs observed should exhibit trajectories
which are relatively unperturbed by galactic and intergalactic magnetic fields. Hence, it is natural to
wonder whether anisotropy begins to emerge at these high energies. Furthermore, if the observed flux

10

850

Energy (eV)

Fig. 6: The compilation of fluorescence-based measurements of the mean X max of various experiments is com-
pared with predictions from EAS simulations using three event generators of hadronic interactions (SIBYLL,
QGSJET, and DPMJET). The QGSJET predictions on the shower- to- shower fluctuations of the depth of maximum
are indicated by the shaded (cross-hatched) area for proton (iron) primaries. From Ref. [19].

suppression is the GZK effect, there is necessarily some distance, 0(100 Mpc), beyond which cosmic
rays with energies near 10 11 GeV will not be seen. Since the matter density within about 100 Mpc is
not isotropic, this compounds the potential for anisotropy to emerge in the UHECR sample. On the one
hand, if the distribution of arrival directions exhibits a large-scale anisotropy, this could indicate whether
or not certain classes of sources are associated with large-scale structures (such as the Galactic plane
or the Galactic halo). On the other hand, if cosmic rays cluster within a small angular region or show
directional alignment with powerful compact objects, one might be able to associate them with isolated
sources in the sky.

CR air shower detectors which experience stable operation over a period of a year or more can have
a uniform exposure in right ascension, a. A traditional technique to search for large-scale anisotropics is
then to fit the right ascension distribution of events to a sine wave with period 2ii/m (m th harmonic) to
determine the components (x, y) of the Rayleigh vector [43]

2 N 2 N

x=— r y^ j w i cos(ra OLi) , V = U] i sin(m a { ) , (2)

i=i i=i

where the sum runs over the number of N events in the considered energy range, J\f = J2iLi w % * s
the normalization factor, and the weights, W{ = are the reciprocal of the relative exposure,

cj, given as a function of the declination, Si [44]. The m th harmonic amplitude of N measurements of
ai is given by the Rayleigh vector length 1Z = [x 2 + y 2 ) 1 / 2 , and the phase is (/? = arctan(y/x).
The expected length of such a vector for values randomly sampled from a uniform phase distribution is
T^o = 2/ \[M . The chance probability of obtaining an amplitude with length larger than that measured is
p(> H) = e~ k °, where ko = 1Z 2 /1Zq. To give a specific example, a vector of length ko > 6.6 would

11

be required to claim an observation whose probability of arising from random fluctuation was 0.0013
(a "3cr" result) [45]. For a given CL, upper limits on the amplitude can be derived using a distribution
drawn from a population characterized by an anisotropy of unknown amplitude and phase

CL= J-

[2 1 ds (Tls\ ( s 2 + K 2 /2 \

where To is the modified Bessel function of the first kind with order zero and a — \/2/J\[ [43].

The first harmonic amplitude of the CR right ascension distribution can be directly related to the
amplitude a^ of a dipolar distribution of the form J(a, S) = Jo(l + add-u), where u and d respectively
denote the unit vector in the direction of an arrival direction and in the direction of the dipole. Setting
m = 1, we can rewrite x, y and J\f as:

x = — \ dS da cos S J(a, 5) u(S) cos a,

2

u min

pSmax P^TT

/ dS da cos S J (a, S) u(S) sin a, (4)

Af = d5 da cos S J (a, S) oj(S) .

In (4) we have neglected the small dependence on right ascension in the exposure. Next, we write the
angular dependence in J (a, 5) as d • u — cos 5 cos (5o cos(o: — ao) + sin 5 sin So, where a^ and 5$ are
the right ascension and declination of the apparent origin of the dipole, respectively. Performing the a
integration in (4) it follows that

Ad±

11 =

where

1 + Bd\\

(5)

J d5cj(5) cos 2 5 f d5 (jj(5) cos 5 sin 5

J d8 uj(S) cos S ' J dS uj(8) cos 6

d|| = a^sin^o is the component of the dipole along the Earth rotation axis, and d± = a^ cos So is
the component in the equatorial plane [46]. The coefficients A and B can be estimated from the data
as the mean values of the cosine and the sine of the event declinations. For example, for the Auger
data sample we have A = (cos S) ~ 0.78 and B = (sin S) ~ —0.45. For a dipole amplitude a^ the
measured amplitude of the first harmonic in right ascension 1Z thus depends on the region of the sky
observed, which is essentially a function of the latitude of the observatory 4ite> and the range of zenith
angles considered. In the case of a small Bd\\ factor, the dipole component in the equatorial plane d± is
obtained as d± ~1Z/ (cos S). The phase tp corresponds to the right ascension of the dipole direction a$.
For a fixed number of arrival directions, the RMS error in the amplitude, Aa^ « 1.57V -1 / 2 , has little
dependence on the amplitude [44]. 1

In Fig. 7 we show upper limits and measurements of d± from various experiments together with
some predictions from UHECR models of both galactic and extragalactic origin. The AGAS A Collabo-
ration reported a correlation of the CR arrival directions to the Galactic plane at the 4<r level [49]. The
energy bin width which gives the maximum fco-value corresponds to the region 10 8 9 GeV - 10 9 3 GeV

AGASA r

where ko = 11.1, yielding a chance probability of p(> 1Z E _ EeV ) ~ 1.5 x 10 . The recent results

1 A point worth noting at this juncture: A pure dipole distribution is not possible because the cosmic ray intensity cannot
be negative in half of the sky. A "pure dipole deviation from isotropy" means a superposition of monopole and dipole, with
the intensity everywhere > 0. An approximate dipole deviation from isotropy could be caused by a single strong source if
magnetic diffusion or dispersion distribute the arrival directions over much of the sky. However, a single source would produce
higher-order moments as well. An example is given in Sec. 2.4.

12

reported by the Pierre Auger Collaboration are inconsistent with those reproted by the AGASA Col-
laboration [48]. If the galactic/extragalactic transition occurs at the ankle, UHECRs at 10 9 GeV are
predominantly of galactic origin and their escape from the Galaxy by diffusion and drift motions are
expected to induce a modulation in this energy range. These predictions depend on the assumed galactic
magnetic field model as well as on the source distribution and the composition of the UHECRs. Two
alternative models are displayed in Fig. 7, corresponding to different geometries of the halo magnetic
fields [50]. The bounds reported by the Pierre Auger Collaboration already exclude the particular model
with an antisymmetric halo magnetic field (A) and are starting to become sensitive to the predictions of
the model with a symmetric field (S). The predictions shown in Fig. 7 are based on the assumption of
predominantly heavy composition in the galactic component [51]. Scenarios in which galactic protons
dominate at 10 9 GeV would typically predict a larger anisotropy. Alternatively, if the structure of the
magnetic fields in the halo is such that the turbulent component dominates over the regular one, purely
diffusion motions may confine light elements of galactic origin up to ~ 10 9 GeV, and may induce an
ankle-like feature at higher energy due to the longer confinement of heavier elements [52]. Typical sig-
natures of such a scenario in terms of large scale anisotropics are also shown in Fig. 7 (dotted line). The
corresponding amplitudes are challenged by the current sensitivity of Auger. On the other hand, if the
transition is taking place at lower energies, say around the second knee, UHECRs above 10 9 GeV are
dominantly of extragalactic origin and their large scale distribution could be influenced by the relative
motion of the observer with respect to the frame of the sources. If the frame in which the UHECR
distribution is isotropic coincides with the CMB rest frame, a small anisotropy is expected due to the
Compton-Getting effect [53]. Neglecting the effects of the galactic magnetic field, this anisotropy would
be a dipolar pattern pointing in the direction ao ~ 168° with an amplitude of about 0.6% [54], close to
the upper limits set by the Pierre Auger Collaboration. The statistics required to detect an amplitude of
0.6% at 99% CL is ~ 3 times the published Auger sample [48].

The right harmonic analyses are completely blind to intensity variations which depend only on
declination. Combining anisotropy searches in a over a range of declinations could dilute the results,
since significant but out of phase Rayleigh vectors from different declination bands can cancel each
other out. An unambiguous interpretation of anisotropy data requires two ingredients: exposure to the
full celestial sphere and analysis in terms of both celestial coordinates [44].

One way to increase the chance of success in finding out the sources of UHECRs is to check
for correlations between CR arrival directions and known candidate astrophysical objects. To calculate
a meaningful statistical significance in such an analysis, it is important to define the search procedure
a priori in order to ensure it is not inadvertently devised especially to suit the particular data set after
having studied it. With the aim of avoiding accidental bias on the number of trials performed in selecting
the cuts, the Auger anisotropy analysis scheme followed a pre-defined process. First an exploratory data
sample was employed for comparison with various source catalogs and for tests of various cut choices.
The results of this exploratory period were then used to design prescriptions to be applied to subsequently
gathered data.

Based on the results of scanning an exposure of 4, 390 km 2 sr yr, a prescription was designed to test
the correlation of events having energies E > 5.5 x 10 10 GeV with objects in the Veron-Cetty & Veron
(VCV) catalog of Active Galactic Nuclei (AGNs). The prescription called for a search of 3.1° windows
around catalog objects with redshifts z < 0.0018. The significance threshold set in the prescription was
met in 2007, when the exposure more than doubled and the total number of events reached 27, with 9
of the 13 events in the post-prescription sample correlating [55,56]. For a sample of 13 events from
an isotropic distribution, the probability that 9 or more correlate by chance with an object in the AGN
catalogue (subject to cuts on the exposure weighted fraction of the sky within the opening angles and the
redshift) is less than 1%. This corresponds to roughly a 2.5<r effect. In the summer of 2008, the HiRes
Collaboration applied the Auger prescription to their data set and found no significant correlation [57].
One has to exercise caution when comparing results of different experiments with potentially different

13

10

i-4

10 19 10 :
Energy [eV]

10

14

10

,15

10

16

10

17

10

18

20

Fig. 7: Upper limits on the anisotropy amplitude of first harmonic as a function of energy from Auger, EAS-TOP,
AGASA, KASCADE and KASCADE-Grande experiments. Also shown are the predictions up to 1 EeV from two
different galactic magnetic field models with different symmetries (A and S), the predictions for a purely galactic
origin of UHECRs up to a few times 10 10 GeV (Gal), and the expectations from the Compton-Getting effect
for an extragalactic component isotropic in the CMB rest frame (C-G XGal). Below 1 EeV, due to variations
of the event counting rate arising from atmospheric effects, the Pierre Auger Collaboration adopted the East/West
method [47], which is two times less efficient but doesn't require correction for trigger efficiency. From Ref. [48]

energy scales, since the analysis involves placing an energy cut on a steeply falling spectrum. More
recently, the Pierre Auger Collaboration published an update on the correlation results from an exposure
of 20, 370 km 2 sr yr (collected over 6 yr but equivalent to 2.9 yr of the nominal exposure/yr of the full
Auger), which contains 69 events with E > 5.5 x 10 10 GeV [58]. A skymap showing the locations of all
these events is displayed in Fig. 8. For a physical signal one expects the significance to increase as more
data are gathered. In this study, however, the significance has not increased. A 3a effect is not neccesarily
cause for excitement; of every 100 experiments, you expect about one 3a effect. Traditionally, in particle
physics there is an unwritten 5a rule for "discovery." One should keep in mind though that in the case of
CR physics we do not have the luxury of controlling the luminosity.

A number of other interesting observations are described in [58], including comparisons with other
catalogs as well as a specific search around the region of the nearest active galaxy, Centaurus A (Cen A).
It is important to keep in mind that these are all a posteriori studies, so one cannot use them to determine
a confidence level for anisotropy as the number of trials is unknown. A compelling concentration of
events in the region around the direction of Cen A has been observed. As one can see in Fig. 9, the
maximum departure from isotropy occurs for a ring of 18° around the object, in which 13 events are
observed compared to an expectation of 3.2 from isotropy. There are no events coming from less than
18° around M87, which is almost 5 times more distant than Cen A and lies at the core of the Virgo
cluster. As shown in Fig. 8 the Auger exposure is 3 times smaller for M87 than for Cen A. Using these
two rough numbers and assuming equal luminosity, one expects 75 times fewer events from M87 than
from Cen A. Hence, the lack of events in this region is not completely unexpected.

14

Fig. 8: Correlation of the arrival directions of UHECR with AGNs from the VCV catalog. The shaded part of the
sky is not visible by Auger. The gray squares are the AGNs within z less than 0.018. The Auger events are shown
with circles. The first 27 events are half filled. The 13 HiRes events are shown with black dots. The thin lines
show the six regions of the sky to which Auger has equal exposure. The wide gray line is the supergalactic plane.
FromRef. [5].

The Centaurus cluster lies 45 Mpc behind Cen A. An interesting question then is whether some
of the events in the 18° circle could come from the Centaurus cluster rather than from Cen A. This does
not appear likely because the Centaurus cluster is farther away than the Virgo cluster and for compara-
blel CR luminosities one would expect a small fraction of events coming from Virgo. Furthermore, the
events emitted by Cen A and deflected by magnetic fields could still register as a correlation due to the
overdense AGN population lying behind Cen A, resulting in a spurious signal [59].

Fig. 9: Cumulative number of events with E > 55 EeV as a function of angular distance from the direction of Cen
A. The bands correspond to the 68%, 95%, and 99.7% dispersion expected for an isotropic flux. From Ref. [58].

15

In summary, the inaugural years of data taking at the Pierre Auger Observatory have yielded
a large, high-quality data sample. The enormous area covered by the surface array together with an
excellent fluorescence system and hybrid detection techniques have provided us with large statistics,
good energy resolution, and solid control of systematic uncertainties. Presently, Auger is collecting some
7, 000 km 2 sr yr of exposure each year, and is expected to run for 2 more decades. New detector systems
are being deployed, which will lower the energy detection threshold down to 10 8 GeV. An experimental
radio detection program is also co-located with the observatory and shows promising results. As always,
the development of new analysis techniques is ongoing, and interesting new results can be expected.

2.2 Origin of ultrahigh energy cosmic rays

It is most likely that the bulk of the cosmic radiation is a result of some very general magneto-hydrodynamic
(MHD) phenomenon in space which transfers kinetic or magnetic energy into cosmic ray energy. The
details of the acceleration process and the maximum attainable energy depend on the particular physical
situation under consideration. There are basically two types of mechanism that one might invoke. The
first type assumes the particles are accelerated directly to very high energy (VHE) by an extended electric
field [60]. This idea can be traced back to the early 1930's when Swann [61] pointed out that betatron
acceleration may take place in the increasing magnetic field of a sunspot. These so-called "one-shot"
mechanisms have been worked out in greatest detail, and the electric field in question is now generally
associated with the rapid rotation of small, highly magnetized objects such as neutron stars (pulsars) or
AGNs. Electric field acceleration has the advantage of being fast, but suffers from the circumstance that
the acceleration occurs in astrophysical sites of very high energy density, where new opportunities for
energy loss exist. Moreover, it is usually not obvious how to obtain the observed power law spectrum in a
natural way, and so this kind of mechanism is not widely favored these days. The second type of process
is often referred to as statistical acceleration, because particles gain energy gradually by numerous en-
counters with moving magnetized plasmas. These kinds of models were mostly pioneered by Fermi [62].
In this case the E~ 2 spectrum emerges very convincingly. However, the process of acceleration is slow,
and it is hard to keep the particles confined within the Fermi engine. In this section we first provide a
summary of statistical acceleration based on the simplified version given in Ref. [63]. For a more detailed
and rigorous discussion, the reader is referred to [64]. After reviewing statistical acceleration, we turn to
the issue of the maximum achievable energy within diffuse shock acceleration and explore the viability
of some proposed UHECR sources.

2.2.1 Fermi acceleration at shock waves

In his original analysis, Fermi [62] considered the scattering of CRs on moving magnetized clouds.
The right panel of Fig. 10 shows a sketch of these encounters. Consider a CR entering into a single
cloud with energy E{ and incident angle 0i with the cloud's direction undergoing diffuse scattering on
the irregularities in the magnetic field. After diffusing inside the cloud, the particle's average motion
coincides with that of the gas cloud. The energy gain by the particle, which emerges at an angle Of
with energy Ef, can be obtained by applying Lorentz transformations between the laboratory frame
(unprimed) and the cloud frame (primed). In the rest frame of the moving cloud, the CR particle has a
total initial energy

E'i = Idoud /^cioud cos 0i) , (6)

where r c i ou d and /3 c ioud = Kioud/c are the Lorentz factor and velocity of the cloud in units of the speed
of light, respectively. In the frame of the cloud we expect no change in energy (E[ = E'j), because
all the scatterings inside the cloud are due only to motion in the magnetic field (so-called collisionless
scattering). There is elastic scattering between the CR and the cloud as a whole, which is much more
massive than the CR. Transforming to the laboratory frame we find that the energy of the particle after

16

its encounter with the cloud is

E f = r cloud E) (1 + /3 cloud cos 6 f ) . (7)
The fractional energy change in the laboratory frame is then

AE _ E f - Ei 1 - /3 d oud cos 0i + /3 doud cos Of - /^ loud cos 0; cos Of

Inside the cloud the CR direction becomes randomized and so (cos Of) = 0. The collisionless scattered
particle will gain energy in a head-on collision (0^ > tt/2) and lose energy by tail-end (0{ < tt/2)
scattering. The net increase of its energy is a statistical effect. The average value of cos 0{ depends on the
relative velocity between the cloud and the particle. The probability P per unit solid angle of having a
collision at angle 0{ is proportional to (v — Kioud cos 0$), where v is the CR speed. In the ultrarelativistic
limit, i.e., v ~ c (as seen in the laboratory frame),

dP

— oc (1 - /3 doud cos 0^) , (9)

so

<cos^) = -^p>. (10)

Now, inserting Eq. (10) into Eq.(8), one obtains for /3 c i ou d <^ 1>

3c2loud/3 - 1 - 4 ^ loud . (id

(AE) _ 1 + /3 c 2 loud /3 4

— — -L r *-~ / —

E 1 - ^loud 3

Note that (AE)/E oc /3 doud , so even though the average magnetic field may vanish, there can still
be a net transfer of the macroscopic kinetic energy from the moving cloud to the particle. However,
the average energy gain is very small, because /3 do ud ^ acceleration process is very similar

to a thermodynamical system of two gases, which tries to come into thermal equilibrium [66]. Corre-
spondingly, the spectrum of CRs should follow a thermal spectrum which might be in conflict with the
observed power-law.

A more efficient acceleration may occur in the vicinity of plasma shocks occurring in astrophysical
environments [67,68]. Suppose that a strong (nonrelativistic) shock wave propagates through the plasma
as sketched in the left panel of Fig. 10. Then, in the rest frame of the shock the conservation relations
imply that the upstream velocity ^ up (ahead of the shock) is much higher than the downstream velocity
^down (behind the shock). The compression ratio r = ^ U p/^down = ^down/^up can be determined by
requiring continuity of particle number, momentum, and energy across the shock; here n up (nd 0W n) is
the particle density of the upstream (downstream) plasma. For an ideal gas the compression ratio can be
related to the specific heat ratio and the Mach number of the shock. For highly supersonic shocks, r =
4 [64]. Therefore, in the primed frame stationary with respect to the shock, the upstream flow approaches
with speed u up = f3 up c = 4/3c/3 and the downstream flow recedes with speed ^ doW n = AiownC = /3c/ 3.
When measured in the stationary upstream frame, the quantity u = u np — ^ doW n = /3c is the speed of
the shocked fluid and u up = /3 s hock is the speed of the shock. Hence, because of the converging flow
- whichever side of the shock you are on, if you are moving with the plasma, the plasma on the other
side of the shock is approaching you with velocity u - to first order there are only head-on collisions for
particles crossing the shock front. The acceleration process, although stochastic, always leads to a gain in
energy. In order to work out the energy gain per shock crossing, we can visualize magnetic irregularities
on either side of the shock as clouds of magnetized plasma of Fermi's original theory. By considering
the rate at which CRs cross the shock from downstream to upstream, and upstream to downstream, one
finds (cos Oi) = —2/3 and (cos Of) = 2/3. Hence, Eq. (8) can be generalized to

( AE ) ~ip = 4 ^np-^down j (12)

E 3 3 c

17

Note this is first order in = u/c, and is therefore more efficient than Fermi's original mechanism.

An attractive feature of Fermi acceleration is its prediction of a power-law flux of CRs. Consider
a test-particle with momentum p in the rest frame of the upstream fluid (see Fig. 10). The particle's
momentum distribution is isotropic in the fluid rest frame. For pitch angles tt/2 < Oi < tt relative to
the shock velocity vector (see Fig. 10) the particle enters the downstream region and has on average
the relative momentum p[l + 2(/3 up — /3down)/3]. Subsequent diffusion in the downstream region 're-
isotropizes' the particle's momentum distribution in the fluid rest frame. As the particle diffuses back
into the upstream region (for pitch angles < Of < tt/2) it has gained an average momentum of
(Ap) jp ~ 4(/3 up — /3down)/3. This means that the momentum gain of a particle per time is proportional
to its momentum,

P = P Again • (13)

On the other hand, the loss of particles from the acceleration region is proportional to their number,

N = -N/t loss . (14)
Therefore, taking the ratio (13)/(14) we first obtain

dN/dp= -aN/p, (15)

and after integration N(p) oc p~ a , with a = t ga { n /t\ oss . If the acceleration cycle across the shock
takes the time At we have already identified At/t ga i n = {Ap) /p ~ 4(/3 up — /?down)/3. For the loss
of relativistic particles one finds At/t\ oss ~ 4/3 up . Therefore, a ~ 3/3 up /(/3 up — /3d 0W n) = 3r/(r — 1),
yielding a ~ 4 for highly supersonic shocks and a > 4 otherwise. The energy spectrum N(E) oc
is related to the momentum spectrum by dEN(E) = Aiip 2 &pN(p) and hence 7 = a — 2 > 2. The
steeply falling spectrum of CRs with 7^3 seems to disfavor supersonic plasma shocks. However, for
the comparison of these injection spectra with the flux of CRs observed on Earth, one has to consider
particle interactions in the source and in the interstellar medium. This can have a great impact on the
shape, as we will discuss in Sees. 2.4 and 2.5.2.

In general, the maximum attainable energy of Fermi's mechanism is determined by the time scale
over which particles are able to interact with the plasma. For the efficiency of a "cosmic cyclotron"
particles have to be confined in the accelerator by its magnetic field B over a sufficiently long time scale

18

compared to the characteristic cycle time. The Larmor radius of a particle with charge Ze increases with
its energy E according to

1 E

Ana ZeB

-¥Gi*)(£)"V

The particle's energy is limited as its Larmor radius approaches the characteristic radial size i? S ource of
the source

E -^ z (^)(^) xio9Gev - <iv '

This limitation in energy is conveniently visualized by the 'Hillas plot' [60] shown in Fig. 1 1 where the
characteristic magnetic field B of candidate cosmic accelerators is plotted against their characteristic size
R. It is important to stress that in some cases the acceleration region itself only exists for a limited period
of time; for example, supernovae shock waves dissipate after about 10 4 yr. In such a case, Eq. (17) would
have to be modifed accordingly. Otherwise, if the plasma disturbances persist for much longer periods,
the maximum energy may be limited by an increased likelihood of escape from the region. A look at
Fig. 11 reveals that the number of sources for the extremely high energy CRs around 10 12 GeV is very
sparse. For protons, only radio galaxy lobes and clusters of galaxies seem to be plausible candidates.
For nuclei, terminal shocks of galactic superwinds originating in the metally-rich starburst galaxies are
potential sources [69]. Exceptions may occur for sources which move relativistically in the host-galaxy
frame, in particular jets from AGNs and gamma-ray bursts (GRBs). In this case the maximal energy
might be increased due to a Doppler boost by a factor ^ 30 or ^ 1000, respectively.

For an extensive discussion on the potential CR-emitting- sources shown in Fig. 11, see e.g. [70].
Two of the most attractive examples are discussed next.

2.2.2 AGNs

AGNs are composed of an accretion disk around a central super-massive black hole and are sometimes
associated with jets terminating in lobes which can be detected in radio. One can classify these objects
into two categories: radio-quiet AGN with no prominent radio emission or jets and radio-loud objects
presenting jets.

Fanaroff-Riley II (FRII) galaxies [71] are the largest known dissipative objects (non-thermal sources)
in the cosmos. Localized regions of intense synchrotron emission, known as "hot spots," are observed
within their lobes. These regions are presumably produced when the bulk kinetic energy of the jets
ejected by a central active nucleus is reconverted into relativistic particles and turbulent fields at a "work-
ing surface" in the head of the jets [72]. Specifically, the speed ^head ~ ^jet [1 + (^e/^jet) 1 ^ 2 ] -1 , with
which the head of a jet advances into the intergalactic medium of particle density n e can be obtained
by balancing the momentum flux in the jet against the momentum flux of the surrounding medium;
where nj e t and Uj e t are the particle density and the velocity of the jet flow, respectively (for relativistic
corrections, see [73]). For n e > nj e t, ^jet > ^head so that that the jet decelerates. The result is the
formation of a strong collisionless shock, which is responsible for particle reacceleration and magnetic
field amplification [74]. The acceleration of particles up to ultrarelativistic energies in the hot spots is the
result of repeated scattering back and forth across the shock front, similar to that discussed in Sec. 2.2.1.
The particle deflection in this mechanism is dominated by the turbulent magnetic field with wavelength
k equal to the Larmor radius of the particle concerned [75]. A self-consistent (although possibly not
unique) specification of the turbulence is to assume that the energy density per unit of wave number of
MHD turbulence is of the Kolmogorov type, I(k) oc A: -5 / 3 , just as for hydrodynamical turbulence [76].
With this in mind, to order of magnitude accuracy using effective quantities averaged over upstream (jet)

19

GT -

MT -

0)

o

• r— I

bJO

a

kT -

T - LHC

G -

mG \-
nG

— I T

neutron
\ stars

i

proton

v,,, — \ \rhite
% *i'v \dwarf

I

,1V\\AGN

'*4

sun GRB\blazar
spot \^

4\ \v

micro
quasar

o-# ^ \ V N

interplanet. SNR
medium

%,galaxy
cipster

= intergal.

g ^l£k lc hal °# 1 medium

starbursfe.
\ wind*.

■ \ \

km Mm GmAU pc kpc Mpc cH 1

characteristic size R

Fig. 11: The "Hillas plot" for various CR source candidates (blue areas). Also shown are jet-frame parameters
for blazers, gamma-ray bursts, and microquasars (purple areas). The corresponding point for the LHC beam is
also shown. The red dashed lines show the lower limit for accelerators of protons at the CR knee (~ 10 6 5 GeV),
CR ankle (~ 10 9,5 GeV) and the GZK suppression (~ 10 10,6 GeV). The dotted gray line is the upper limit
from synchrotron losses and proton interactions in the cosmic photon background (R ^> 1 Mpc). The grey area
corresponds to astrophysical environments with extremely large magnetic field energy that would be gravitationally
unstable. FromRef. [65].

and downstream (hot spot) conditions (considering that downstream counts a fraction of 4/5 [75]) the
acceleration timescale at a shock front is found to be [77]

AGN 20Z)|| (E)

u

where

DAE)

2c

jet

eB J

(18)

(19)

is the Kolmogorov diffusion coefficient, U is the ratio of turbulent to ambient magnetic energy density
in the region of the shock (of radius R), and B is the total magnetic field strength.

20

The subtleties surrounding the conversion of particles kinetic energy into radiation provide ample
material for discussion. The most popular mechanism to date relates 7-ray emission to the development
of electromagnetic cascades triggered by secondary photomeson products that cool instanstaneously via
synchrotron radiation. The characteristic single photon energy in synchrotron radiation emitted by an
electron is

^ n = (£) 1/2 £^f§ ~ 5 - 4 x 10-2 ^° ^/ Eey ) 2 Tev • ™

For a proton this number is (m p /m e ) 3 ~6x 10 9 times smaller.

The acceleration process will then be efficient as long as the energy losses by synchrotron radiation
and/or photon-proton interactions do not become dominant. The synchrotron loss time for protons is
given by [78]

Ts y n ~ 2TD2 ' ( 21 )

where or and T = E/(m p c 2 ) are the Thomson cross section and Lorentz factor, respectively. Con-
sidering an average cross section a 1P [79] for the three dominant pion-producing interactions, jp —>
p7T°, 7p —> n7r + , 7p —> p7T + 7T~ , the time scale of the energy losses, including synchrotron and photon
interaction losses, reads [80]

6tt mi C 3 - T SV n

E' 1 = syn A , (22)

loss ~ a T m 2 B 2 (1 + Aa) 1 + Aa

where a stands for the ratio of photon to magnetic energy densities and A gives a measure of the relative
strength of ^p interactions versus the synchrotron emission. Note that the second channel involves the
creation of ultrarelativistic neutrons that can readily escape the system. For typical hot spot conditions,
the number density of photons per unit energy interval follows a power-law spectrum

AGN/ a J (N /CJ ) (CJ/CJO)- 2 U Q <UJ<UJ*

^ W-\ otherwise ^

where No is the normalization constant and uo and uo* correspond to radio and gamma rays energies,
respectively. The ratio of photon to magnetic energy density is then

a ~ B 2 /8tt (M)
and A is only weakly dependent on the properties of the source

A= a JL (m p /m e ) 2 ^^n> lQx 1q5 ^ 20Q (25)
ax m(o;*/a;o)

The maximum attainable energy can be obtained by balancing the energy gains and losses [81]

E 20 = 1.4 x 10 5 B-W f% (1 + Aa)^ , (26)

where E = 10 20 E20 eV and R = R\^ pc 1 kpc. It is of interest to apply the acceleration conditions to the
nearest AGN.

At only 3.4 Mpc distance, Cen A is a complex FRI radio-loud source identified at optical frequen-
cies with the galaxy NGC 5128 [82]. Radio observations at different wavelengths have revealed a rather
complex morphology shown in Fig. 12. It comprises a compact core, a jet (with subluminal proper mo-
tions /3jet ~ 0.5 [85]) also visible at X-ray frequencies, a weak counterjet, two inner lobes, a kpc-scale
middle lobe, and two giant outer lobes. The jet would be responsible for the formation of the northern
inner and middle lobes when interacting with the interstellar and intergalactic media, respectively. There

21

Fig. 12: Optical image of Cen A (UK 48-inch Schmidt) overlaid with radio contours (black, VLA [83]). H.E.S.S.
VHE best fit position with la statistical errors (blue cross) and VHE extension 95% CL upper limit (white dashed
circle) are also shown. From Ref. [84].

appears to be a compact structure in the northern lobe, at the extrapolated end of the jet. This structure
resembles the hot spots such as those existing at the extremities of FRII galaxies. However, at Cen A it
lies at the side of the lobe rather than at the most distant northern edge, and the brightness contrast (hot
spot to lobe) is not as extreme [86]. Estimates of the radio spectral index of synchrotron emission in
the hot spot and the observed degree of linear polarization in the same region suggests that the ratio of
turbulent to ambient magnetic energy density in the region of the shock is U ~ 0.4 [87]. The broadband
radio- to-X-ray jet emission yields an equipartition magnetic field B^q ~ 100 [88]. 2 The radio- visible
size of the hot spot can be directly measured from the large scale map R\^ pc ~ 2 [89]. The actual size
can be larger because of uncertainties in the angular projection of this region along the line of sight. 3 Re-
placing these fiducial values in (17) and (26) we conclude that if the ratio of photon to magnetic energy
density a < 0.4, it is plausible that Cen A can accelerate protons up E « 2 x 10 11 GeV.

EGRET observations of the gamma ray flux for energies > 100 MeV allow an estimate L 7 ~
10 41 erg s _1 for Cen A [90]. This value of L 7 is consistent with an earlier observation of photons in the
TeV-range during a period of elevated X-ray activity [91], and is considerably smaller than the estimated
bolometric luminosity Lbol ~ 10 43 erg s _1 [82]. Recent data from H.E.S.S. have confirmed Cen A as a
TeV 7-ray emitting source [84]. Extrapolating the spectrum measured with EGRET in the GeV regime

2 The usual way to estimate the magnetic field strength in a radio source is to minimize its total energy. The condition of
minimum energy is obtained when the contributions of the magnetic field and the relativistic particles are approximately equal
(equipartition condition). The corresponding B-field is commonly referred to as the equipartion magnetic field.

3 For example, an explanation of the apparent absence of a counterjet in Cen A via relativistic beaming suggests that the
angle of the visible jet axis with respect to the line of sight is at most 36° [86], which could lead to a doubling of the hot spot
radius. It should be remarked that for a distance of 3.4 Mpc, the extent of the entire source has a reasonable size even with this
small angle.

22

to VHEs roughly matches the H.E.S.S. spectrum, though the softer end of the error range on the EGRET
spectral index is preferred. More recent data from Fermi-LAT established that a large fraction (> 1/2)
of the total > 100 MeV emission from Cen A emanates from the lobes [92]. For values of B in the
100 fiG range, substantial proton synchrotron cooling is suppressed, allowing production of high energy
electrons through photomeson processes. The average energy of synchrotron photons scales as [93]

{E^ n ) -0.29£; yn , (27)

and therefore, to account for the observed TeV photons Cen A should harbor a population of ultra-
relativistic electrons, with E e ~ 10 9 GeV. We further note that this would require the presence of protons
with energies between one and two orders of magnitude larger, since the electrons are produced as
secondaries.

2.2.3 GRBs

GRBs are flashes of high energy radiation that can be brighter, during their brief existence, than any other
source in the sky. The bursts present an amazing variety of temporal profiles, spectra, and timescales [94].
Our insights into this phenomenon have been increased dramatically by BATSE observations of over
2000 GRBs, and more recently, by data from SWIFT.

There are several classes of bursts, from single-peaked events, including the fast rise and expo-
nential decaying (FREDs) and their inverse (anti-FREDs), to chaotic structures [95]. There are well
separated episodes of emission, as well as bursts with extremely complex profiles. Most of the bursts are
time asymmetric, but some are symmetric. Burst timescales range from about 30 ms to several minutes.

The GRB angular distribution appears to be isotropic, suggesting a cosmological origin [96]. Fur-
thermore, the detection of "afterglows" — delayed low energy (radio to X-ray) emission — from GRBs
has confirmed this via the redshift determination of several GRB host-galaxies [97].

The 7-ray luminosity implied by cosmological distances is astonishing: L 1 ~ 10 52 erg/s. The
most popular interpretation of the GRB -phenomenology is that the observable effects are due to the
dissipation of the kinetic energy of a relativistic expanding plasma wind, a "fireball" [98]. Although the
primal cause of these events is not fully understood, it is generally believed to be associated with the
core collapse of massive stars (in the case of long duration GRBs) and stellar collapse induced through
accretion or a merger (short duration GRBs) [99].

The very short timescale observed in the light curves indicates an extreme compactness (i.e. dis-
tance scale comparable to a light-ms: ro ~ 10 7 cm) that implies a source which is initially opaque
(because of 77 pair creation) to 7-rays

r 77 ~ r n« RB a T ~ ~ 10 15 , (28)

1 1 j 47T ro c e 7

where n^ RB is number density of photons at the source and e 7 ~ 1 MeV is the characteristic photon
energy.

The high optical depth creates the fireball: a thermal plasma of photons, electrons, and positrons.
The radiation pressure on the optically thick source drives relativistic expansion (over a time scale ro / c),
converting internal energy into the kinetic energy of the inflating shell. As the source expands, the optical
depth is reduced. If the source expands with a Lorentz factor T, the energy of photons in the source frame
is smaller by a factor T compared to that in the observer frame, and most photons may therefore be below
the pair production threshold. Baryonic pollution in this expanding flow can trap the radiation until most
of the initial energy has gone into bulk motion with Lorentz factors of T ~ 10 2 — 10 3 [100]. The
kinetic energy can be partially converted into heat when the shell collides with the interstellar medium
or when shocks within the expanding source collide with one another. The randomized energy can then
be radiated by synchrotron radiation and inverse Compton scattering yielding non-thermal bursts with

23

timescales of seconds at large radii, r > 10 12 cm, beyond the Thompson sphere. Charged particles may
be efficiently accelerated to ultrahigh energies in the fireball's internal shocks, hence GRBs are often
considered as potential sources of UHECRs [101].

Coburn and Boggs [102] reported the detection of polarization, a particular orientation of the
electric-field vector, in the 7-rays observed from a burst. The radiation released through synchrotron
emission is highly polarized, unlike in other previously suggested mechanisms such as thermal emission
or energy loss by relativistic electrons in intense radiation fields. Thus, polarization in the 7-rays from a
burst provides direct evidence in support of synchrotron emission as the mechanism of 7-ray production
(see also [103]).

Following Hillas criterion, the Larmor radius tl should be smaller than the largest scale Zgrb
over which the magnetic field fluctuates, since otherwise Fermi acceleration will not be efficient. One
may estimate Zgrb as follows. The comoving time, i.e. the time measured in the fireball rest frame, is
t = r/Tc. Hence, the plasma wind properties fluctuate over comoving scale length up to /grb ~ r/T,
because regions separated by a comoving distance larger than r/T are causally disconnected. Moreover,
the internal energy is decreasing because of the expansion and thus it is available for proton acceleration
(as well as for 7-ray production) only over a comoving time t. The typical acceleration time scale is
then [101]

^ B ~ T ± • (29)

Equation (29) sets a constraint on the required comoving magnetic field strength, and the Larmor radius
tl — E' /eB — E/FeB, where E' — E/Y is the proton energy measured in the fireball frame. This
constraint sets a lower limit to the magnetic field carried by the wind, which may be expressed as

<B > Q M T h^k^ (3Q)
Ce £52

where T = 10 2 5 r2.5, L 1 = 10 52 Ls2 erg s _1 . Here, (b is the fraction of the wind energy density which
is carried by the magnetic field, Anr 2 T 2 (B 2 / '8tt) = (bL, and ( e is the fraction of wind energy carried
by shock accelerated electrons. Note that because the electron energy is lost radiatively, L 1 « ( e L.

The dominant energy loss process in this case is synchrotron cooling. Therefore, the condition
that the synchrotron loss time of Eq. (21) be smaller than the acceleration time sets the upper limit on the
magnetic field strength [101]

" ° w m T 2.5 ^20

Since the equipartition field is inversely proportional to the radius r, this condition may be satisfied
simultaneously with (30) provided that the dissipation radius is large enough, i.e.

r > 10 12 T^ 2 E% cm. (32)

The high energy protons also lose energy in interaction with the wind photons (mainly through pion
production). It can be shown, however, that this energy loss is less important than the synchrotron
energy loss [101].

In summary, a dissipative ultra-relativistic wind, with luminosity and variability time implied by
GRB observations, satisfies the constraints necessary to accelerate protons to energy > 10 11 GeV, pro-
vided that T > 100, and the magnetic field is close to equipartition with electrons.

B < 3 x 10 5 r| 5 E~ 2 G. (31)

2.3 Energy losses of baryonic cosmic rays on the pervasive radiation fields
2.3.1 Opacity of the CMB to UHECR protons

Ultrahigh energy protons degrade their energy through Bethe-Heitler (BH) pair production (pj —>
pe + e~) and photopion production (pj — >► ttN), each successively dominating as the proton energy in-
creases. The fractional energy loss due to interactions with the cosmic background radiation at a redshift

24

z = is determined by the integral of the nucleon energy loss per collision multiplied by the probability
per unit time for a nucleon collision in an isotropic gas of photons [104]. This integral can be explicitly
written as follows,

IdE c ^ [ u ™ 1 . . nJu)

~Elt = 2^^ J * ^ i r/2T du ~W ' (33)

where uo r is the photon energy in the rest frame of the nucleon, and yj is the inelasticity, i.e. the average
fraction of the energy lost by the photon to the nucleon in the laboratory frame for the jth reaction
channel. (Here the laboratory frame is the one in which the CMB is at a temperature ^ 2.7 K.) The sum
is carried out over all channels, n 7 (uj)du stands for the number density of photons with energy between
uj and duj, aj(uj r ) is the total cross section of the jth interaction channel, T is the usual Lorentz factor of
the nucleon, and uj m is the maximum energy of the photon in the photon gas.

Pair production and photopion production processes are only of importance for interactions with
the 2.7 K blackbody background radiation [105]. Collisions with optical and infrared photons give a
negligible contribution. Therefore, for interactions with a blackbody field of temperature T, the photon
density is that of a Planck spectrum, so the fractional energy loss is given by

where uo j is the threshold energy for the jth reaction in the rest frame of the nucleon.

At energies E <C m e m p /kT = 2.1 x 10 9 GeV (i.e., uj r /m e - 2 C 1), when the reaction takes
place on the photons from the high energy tail of the Planck distribution, the fraction of energy lost in
one collision and the cross section can be approximated by the threshold values

V m = 2 ^ , (35)

I Up

and

W e ~ 2 ) ' (36)

where a is the fine structure constant and ro is the classical radius of the electron [105]. The fractional
energy loss due to pair production for E < 10 9 GeV is then,

1 (dE\ lQcm e 2 ( kT\ 3 (FkT\ 2 / m e \

~e \ii) Bn = — m P ar ° \m) UtJ ex n-fkf)- (37)

At higher energies (E > 10 10 GeV) the characteristic time for the energy loss due to pair production
is t ~ 5 x 10 9 yr [106]. In this energy regime, the photopion reactions —> pn and pj —> n+n on
the tail of the Planck distribution give the main contribution to proton energy loss. The cross sections of
these reactions are well known and the kinematics is simple.

Photopion production turns on at a photon energy in the proton rest frame of 145 MeV with a
strongly increasing cross section at the A (1232) resonance, which decays into the one pion channels
7T + n and 7T°p. With increasing energy, heavier baryon resonances occur and the proton might reappear
only after successive decays of resonances. The most important channel of this kind is p^ A ++ 7r~
with intermediate A ++ states leading finally to A ++ — >► p?r + . A ++ examples in this category are the
A(1620) and A(1700) resonances. The cross section in this region can be described by either a sum or
a product of Breit-Wigner distributions over the main resonances produced in 7 collisions considering
tvN, tvtvN and KA (A Ntv) final states [107]. At high energies, 3.0 GeV < u r < 183 GeV, the

3

25

E\ = vi^A_ ^ (40)

CERN-HERA and COMPAS Groups have made a fit to the p^ cross section [108]. The parameterization
is

a n (u r ) = A + B In 2 (^) + C In (^) mb , (38)

where A = 0.147 ± 0.001, B = 0.0022 ± 0.0001, and C = -0.0170 ± 0.0007. In this energy range,
the cr to tai(^7) is to a good approximation identical to cr to tai(P7)-

We turn now to the kinematics of photon-nucleon interactions. The inelasticity y n depends not
only on the outgoing particles but also on the kinematics of the final state. Nevertheless, averaging over
final state kinematics leads to a good approximation of y n . The cm. system quantities (denoted by
*) are determined from the relativistic invariance of the square of the total 4-momentum p^ of the
photon-proton system. This invariance leads to the relation

s = (cj* + E*) 2 = m 2 p + 2m p u r . (39)

The cm. system energies of the particles are uniquely determined by conservation of energy and mo-
mentum. For reactions mediated by resonances one can assume a decay, which in the cm. frame is
symmetric in the forward and backward directions with respect to the collision axis (given by the in-
coming particles). For instance, we consider single pion production via the reaction p^ —> A —> pn.
Here,

(s + m\ - m 2 )
2^~s

Thus, the mean energy of the outgoing proton is

/ F *finah = (s + m 2 A - ml) (m| + m 2 - ml)

[ p ' 2^m A 2m A ' K }

/pfinalv E (g ~ ml + mj) {m\ - m 2 + m 2 )
{E ? ) = 7 2^ 2^ • (42)

The mean inelasticity y n = 1 — ( (i? final ) /E) of a reaction that provides a proton after n resonance decays
can be obtained by straightforward generalization of Eq. (42), and is given by

(43)

where m# denotes the mass of the i tn resonant system of the decay chain, mM the mass of the associ-
ated meson, trr = is the total energy of the reaction in the cm., and rriR n the mass of the nucleon.
For multi-pion production the case is much more complicated because of the non-trivial final state kine-
matics. However, it is well established experimentally [109] that, at very high energies (y/s > 3 GeV),
the incoming particles lose only one-half their energy via pion photoproduction independently of the
number of pions produced, y n ~ 1/2. This is the "leading particle effect".

For < 2 GeV, the best maximum likelihood fit to Eq. (34) with the exponential behavior

or in the lab frame

1 fdE
E l¥

Aexp[-B/E], (44)

derived from the values of cross section and fractional energy loss at threshold, gives [110]

A = (3.66 ± 0.08) x 10" 8 yr -1 , B = (2.87 ± 0.03) x 10 11 GeV . (45)

The fractional energy loss due to production of multipion final states at higher cm. energies (y/s >
3 GeV) is roughly a constant,

) =C= (2.42 ± 0.03) x 10~ 8 yr" 1 . (46)
E \dt J n

26

1C?° 1C? 2

E [eV]

Fig. 13: Energy attenuation length of protons in the intergalactic medium. Note that after a distance of ~ 100 Mpc,
or propagation time ~ 3 x 10 8 yr, the mean energy is essentially independent of the initial energy of the protons,
with a critical energy around 10 11 GeV. From Ref. [110].

From the values determined for the fractional energy loss, it is straightforward to compute the energy
degradation of UHECRs in terms of their flight time. This is given by,

At - Ei (B/E) + Ei (B/Eq) = , for 10 10 GeV < E < 10 12 GeV , (47)

and

E(t) = E exp[- C t ] , for E > 10 12 GeV , (48)

where Ei is the exponential integral [111]. Figure 13 shows the proton energy degradation as a function
of the mean flight distance. Notice that, independent of the initial energy of the nucleon, the mean energy
values approach 10 11 GeV after a distance of « 100 Mpc.

2.3.2 Photonuclear interactions

The relevant mechanisms for the energy loss that extremely high energy nuclei suffer during their trip
to Earth are: Compton interactions, pair production in the field of the nucleus, photodisintegration, and
hadron photoproduction. The Compton interactions have no threshold energy. In the nucleus rest-frame,
pair production has a threshold at ^ 1 MeV, photodisintegration is particularly important at the peak of
the GDR (15 to 25 MeV), and photomeson production has a threshold energy of ~ 145 MeV.

Compton interactions result in only a negligibly small energy loss for the nucleus given by [112]

27

where p 1 is the energy density of the ambient photon field in eV cm -3 , E is the total energy of the
nucleus in eV, and Z and A are the atomic number and weight of the nucleus. The energy loss rate
due to photopair production is Z 2 /A times higher than for a proton of the same Lorentz factor [113],
whereas the energy loss rate due to photomeson production remains roughly the same. The latter is
true because the cross section for photomeson production by nuclei is proportional to the mass number
A [114], while the inelasticity is proportional to 1/A However, it is photodisintegration rather than
photopair and photomeson production that determines the energetics of UHECR nuclei. During this
process some fragments of the nuclei are released, mostly single neutrons and protons. Experimental
data of photonuclear interactions are consistent with a two-step process: photoabsorption by the nucleus
to form a compound state, followed by a statistical decay process involving the emission of one or more
nucleons.

Following the conventions of Eq. (33), the disintegration rate with production of i nucleons is
given by [115]

1 f°° n (uj) f 2Vuj

R Ai=^2 J du 1 J1 J d0J r U r (7 A i{Ur) (50)

with cjAi the cross section for the interaction.

The photoabsorption cross section roughly obeys the Thomas-Reiche-Kuhn (TRK) dipole sum

rule

f°° N Z

E d = / a{uj r )duj r = 59.8 —— MeVmb , (51)
Jo A

where N = A — Z is the number of neutrons. (Indeed, this integral is experimentally ^ 20 — 30% larger,

e.g. for 56 Fe, 1, 020 mb-MeV for the left hand side, 22% larger than the right hand side [116].) These

cross sections contain essentially two regimes. At uo r < 30 MeV there is the domain of the GDR where

disintegration proceeds mainly by the emission of one or two nucleons. A Gaussian distribution in this

energy range is found to adequately fit the cross section data [112]. At higher energies, the cross section

is dominated by multinucleon emission and is approximately flat up to u r ~ 150 MeV. Specifically,

- 2) O(30 - ay) e -2K-eo0 2 /A? /iSd e(a; r - 30)

° M = WA~i + 120 ' (52)

for i=\, 2, and

/iS d e(a; r -30) ....

°M = ' (53)

for i > 2 [112]. Here, W is a normalization factor given by

W = g) 2 [$(^2(30 - e 0l )/A z ) + $(^2(6o, - 2)/A<)] ,

<3>(x) is the error function, and Q(x) the Heaviside step function. The dependence of the width A^, the
peak energy eoi, the branching ratio fi, and the dimensionless integrated cross section ^ are given in
Ref. [112] for isotopes up to 56 Fe.

The photon background relevant for nucleus disintegration consists essentially of photons of the
2.7 K CMB. The background of optical radiation turns out to be of (almost) no relevance for UHECR
propagation. The cosmic infrared background (CIB) radiation [117]

^p- = 1.1 x ID" 4 f^-)" 2 ' 5 cm-SeV" 1 , (54)
duo VeV/

only leads to sizeable effects far below 10 11 GeV and for time-scales O (10 17 s) [118].

By substituting Eqs. (52) and (53) into Eq. (50) the photodisintegration rates on the CMB can be
expressed as integrals of two basic forms. The first one is

*i5/r

2r 2 7T 2 ft 3 C 2

ri-o/L roc
Jl/Y J15/T

(55)

28

where the functions J and J' are given by the expressions,

J =

+

eoiAi \ <S>(V2{2Tuj - eoi)/A<) + <S>(V2(e 0i - 2)/A<

(56)

A.

-2((e *-2)/A,) 2 _ -2((2ru;-e i)/Ai

'}•

and

J' = J-eoiAi [$(^2(30 - eoi)/A<) + *( V2(e 0i - 2)/A<

2

(57)

+

9<

e -2((e ,-2)/A,) 2 _ e -2((30-e 0l )/A,

The second basic integral is of the form

I 2 = (tt^VJ-Va

ii 5 /r e"/ fcT - 1 ~ V ~r) L/r

e uj/kT i

(58)

With this in mind, Eq. (50) can be re- written as [1 19]

^J f kT\n(l-e- 15 / rkT )

-±e-*&M (£) VZ f <S 3 + | JCkT [HI ~ e~^ kT ) - ln(l - e~V™)

TT \ 1/2 A?

A

+

120

r 2 5 4 + 15 2 fcrin(l-e- 15 / rfer )

(59)

with A, Si, and K as given in Table 2. Summing over all the possible channels for a given number of
nucleons, one obtains the effective nucleon loss rate Ra = iR-Ai- The effective nucleon loss rate for
light elements, as well as for those in the carbon, silicon and iron groups can be scaled as [1 12]

dA
~dt

dA
~dt

56 Fe

#56

A
56

with the photodisintegration rate (59) parametrized by [120]

R 56 (T) = 3.25 x 10- 6 r- a643 exp(-2.15 x 10 10 /r) s" 1
forT G [1.0 x 10 9 ,36.8 x 10 9 ], and

i? 5 e(r) = i.59 x io _i2 r-

forT G [3.68 x 10 10 , 10.0 x 10 10 ].

-12-T-0.0698 -1

(60)

(61)

(62)

EXERCISE 2.1 Approximating the cross section in Eq. (52) by the single pole of the Narrow-
Width Approximation [121]

ga{oj') = ir&o r ° ( f R 8(u>' - luq) , (63)

show that for interactions with the CMB photons

Ra

4r 2 vr

In (l - e-<l 2TT )

(64)

29

Table 2: Series and functions of Eq. (59).

A

Si

E7 =1 kTf 1 exp[B z ] + 15V8/A0 - #(£ + >/8/A<)}

s 2

E7 =1 ^i" 1 exp{- J yrfcT}[$( v / 2(2 - eoiJ/Ai) - $(^(30 - eoO/AO]

S3

1 exp[^] + 15v/8/Ai) - $(B + VS/Ai)}

S 4

El, exp{-15j/rfcT}[(fcT/j)(15/r) 2 + (fcT/j) 2 (15/r) + {kT/jf]

B

jAi/FkTVS2 - 2e 0i /V2Ai

K

e 0i Ai <S>(V2(e 0i - 2)/A;) + (A;/2) 2 exp{-2(£ 0i - 2) 2 /A 2 }

where a /A = 1.45 x lCT 27 cm 2 , r GDR = 8 MeV, and e' = 42.65A~ 21 (0.925A 2 433 ) MeV, for
A > 4 (A < 4) [122]. Verify that for 56 Fe this solution agrees to within 20% with the parametrization
given in Eq. (61).

For photodisintegration, the averaged fractional energy loss results equal the fractional loss in
mass number of the nucleus, because the nucleon emission is isotropic in the rest frame of the nucleus.
During the photodisintegration process the Lorentz factor of the nucleus is conserved, unlike the cases
of pair production and photomeson production processes which involve the creation of new particles that
carry off energy. The total fractional energy loss is then

_ 1 dE _ 1 dY R

~E^tt~f^tt + A' ( }

For u r < 145 MeV the reduction in T comes from the nuclear energy loss due to pair production. The 7-
ray momentum absorbed by the nucleus during the formation of the excited compound nuclear state that
precedes nucleon emission is (9(10 -2 ) times the energy loss by nucleon emission [123]. For T > 10 10
the energy loss due to photopair production is negligible, and thus

E{t) - 938 A(t) T MeV

Eq exp

56

(66)

Figure 14 shows the energy of the heaviest surviving nuclear fragment as a function of the prop-
agation time, for initial iron nuclei. The solid curves are obtained using Eq. (66), whereas the dashed
and dotted-dashed curves are obtained by means of Monte Carlo simulations [118]. One can see that
nuclei with Lorentz factors above 10 10 cannot survive for more than 10 Mpc. For these distances, the
approximation given in Eq. (66) always lies in the region which includes 95% of the Monte Carlo simula-
tions. When the nucleus is emitted with a Lorentz factor Tq < 5 x 10 9 , pair production losses start to be
relevant, significantly reducing the value of T as the nucleus propagates distances of (9(100 Mpc). The
effect has a maximum for Tq « 4 x 10 9 but becomes small again for Tq < 10 9 , for which appreciable
effects only appear for cosmological distances (> 1000 Mpc), see for instance [118].

Note that Eq. (66) imposes a strong constraint on the location of nucleus-sources: less than 1%
of iron nuclei (or any surviving fragment of their spallations) can survive more than 3 x 10 14 s with an
energy > 10 1L5 GeV For straight line propagation, this represents a maximum distance of ^ 3 Mpc.

2.4 Diffuse propagation of protons in a magnetized Local Supercluster

In addition to the interactions with the radiation fields permeating the universe, baryonic CRs suffer
deflection and delay in magnetic fields, effects which can camouflage their origins. For example, the
regular component of the Galactic magnetic field can distort the angular images of CR sources: the flux

30

1e+19

1e+13 1e+14 1e+15 1e+16 1e+17

t[s]

Fig. 14: The energy of the surviving fragment (To = 4 x 10 9 , To = 2 x 10 10 ) vs. propagation time obtained
using Eq. (66) is indicated with a solid line. Also included is the energy attenuation length obtained from Monte
Carlo simulations with (dashed) and without (dotted-dashed) pair creation production, for comparison. The region
between the two dotted lines includes 95% of the simulations. This gives a clear idea of the range of values which
can result from fluctuations from the average behaviour. It is important to keep in mind that a light propagation
distance of 1.03 x 10 14 s corresponds to 1 Mpc. From Ref. [124].

may appear dispersed around the source or globally translated in the sky with rather small dispersion, viz.
the deflection for CRs of charge Ze and energy E should not exceed - 10° Z(4x 10 10 GeV/E) [125].

One interesting possibility to explain the observed near-isotropy of arrival directions is to envis-
age a large scale extragalactic magnetic field that can provide sufficient bending to the CR trajectories.
Surprisingly little is actually known about the extragalactic magnetic field strength. There are some
measurements of diffuse radio emission from the bridge area between the Coma and Abell superclus-
ters [126], which under assumptions of equipartition allows an estimate of (9(0.2 — 0.6) /xG for the
magnetic field in this region. Fields of O(fiG) are also indicated in a more extensive study of 16 low
redshift clusters [127]. It is assumed that the observed £> -fields result from the amplification of much
weaker seed fields. However, the nature of the initial week seed fields is largely unknown. There are two
broad classes of models for seed fields: cosmological models, in which the seed fields are produced in
the early universe, and astrophysical models, in which the seed fields are generated by motions of the
plasma in (proto)galaxies. Of particular interest here is the second class of models. If most galaxies lived
through an active phase in their history, magnetized outflows from their jets and winds would efficiently
pollute the extragalactic medium. The resulting £>-field is expected to be randomly oriented within cells
of sizes below the mean separation between galaxies, A# < 1 Mpc.

Extremely weak unamplified extragalactic magnetic fields have escaped detection up to now. Mea-
surements of the Faraday rotation in the linearly polarized radio emission from distant quasars [128]
and/or distortions of the spectrum and polarization properties in the CMB [129, 130] imply upper limits
on the extragalactic magnetic field strength as a function of the reversal scale. It is important to stress that
Faraday rotation measurements (RM) sample extragalactic magnetic fields of any origin (out to quasar
distances), while the CMB analyses set limits only on primordial magnetic fields. The RM bounds de-
pend significantly on assumptions about the electron density profile as a function of the redshift. When
electron densities follow that of the Lyman-a forest, the average magnitude of the magnetic field receives
an upper limit of B ~ 10 -9 G for reversals on the scale of the horizon, and B ~ 10 -8 G for reversal

31

scales on the order of 1 Mpc [131]. As a statistical average over the sky, an all pervading extragalactic
magnetic field is constrained to be [132]

B<3x 1(T 7 (Q^/Om)- 1 (Zi/0.72) (A^/Mpc) 1 / 2 G , (67)

where Q^h 2 ~ 0.02 is the baryon density and h ~ 0.72 is the present day normalized Hubble expansion
rate. This is a conservative bound because has contributions from neutrons and only electrons in
ionized gas are relevant to Faraday rotation.

In the spirit of [133, 134], very recently we proposed that neutron emission from Cen A could
dominate the observed CR flux above the GZK suppression [135]. Neutrons that are able to decay
generate proton diffusion fronts in the intergalactic turbulent magnetic plasma. In our calculations we
assume a strongly turbulent magnetic field: B = 50 nG, Xb ~ 1 Mpc, and largest turbulent eddy
£ ~ 2tt\b [136]. For energies above the GZK supression, Xb < tl < £ and so the diffusion coefficient
is given by the Bohm formula [75]

The evolution of the proton spectrum is driven by the so-called "energy loss-diffusion equation"

dn(E,r,t) _ d[b(E)n(E , r , t)}

(68)

dt dE

+ V[D(E, r, t) Vn(E, r, *)] + Q(E, t) S 3 (x) , (69)

where n(E, r, 0) = No S 3 (x). Here, b(E) = dE/dt is the mean rate at which particles lose energy
and Q(E,t) is the number of protons per unit energy and per unit time generated by the source. For
the situation at hand, D(E,r,t) = D{E) and hence the second term becomes D(E)V 2 n(E,r,t).
Idealizing the emission to be uniform with a rate dNo /dt = N tot / r, we have

Q(E,t) = ^ [Q(t - ton) - 0(t - toff)], (70)
r

where / Q(E, r', t') d 3 x' dt' = N tot , O is the Heaviside step function, and t on (t g) is the time since the
engine turned on (off) its CR production, £ Q ff — ^on = t. For the energy region of interest, the expected
time delay of the diffuse protons, raday ~ d 2 /D(E), is significantly smaller than the characteristic time
scale for photopion production derived in Eq. (44).

If the energy loss term is neglected, the solution to Eq. (69) reads,

n(E, r,t) = j dt' J d 3 x' G(r -r',t- t') Q(E, r', t'), (71)

where

G(r-r',t-t') = [47rD(t-t')]" 3/2 ©(* - O exp{-(r - r') 2 /4D(t - t')} , (72)

is the Green function [137]. The density of protons at the present time t of energy E at a distance r from
Cen A, which is assumed to be continuously emitting at a constant spectral rate dNo / dE dt from time
i on until the present, is found to be [134]

dn(E, r, t) _ dN 1 [ l ^ e -r 2 /4£>(t-f )

dE dE dt [4ttL>(£)] 3 / 2 J tm {t - t'fl 2

r 2 d

)r J vi v 3

dN 1 r dv _ 1/v

dEdt 4n 3 / 2 D(E)r J V1 v 3 ^
dNo 1

dE dt 4irD(E)r

I(x) , (73)

32

where we have used the change of variables

U =4W^)=l (74)

with x = 4DT on /r 2 , T on = t — t on , and

du

1 f°°

I{x) = -f= / e~ u . (75)

For T on — >► oc, the density approaches its time-independent equilibrium value n eq .

As a result of this diffusion the J oc E~ A behavior of the observed CR spectrum reflects a
dNo/dEdt oc E~ 3 injection in the region of the source cutoff. For S = 3000 km 2 detector like Auger,
the neutron rate is

f

dt And? J El dEdt

where \{E) ~ 9.2 x 10" 3 E^eV Mpc is the neutron decay length. For the energy interval between E\ =
55 EeV and E2 = 150 EeV, we calculate the normalization factor using the observation of 2 neutrons
in 3 yr of the nominal exposure/yr of Auger. We then use this normalization factor to calculate the
luminosity of the source in the above energy interval. We find L^ E2) = 0.86 x 10 40 erg/s [135]. Next,
we assume continuity of the spectrum at E\ as it flattens to E~ 2 . Taking the lower bound on the energy
to be Eq = 10 EeV, we can then fix the luminosity for this interval and find L^ ,El ^ = 2.3 x 10 40 erg/s.

Adding these, we find the (quasi) bolometric luminosity to be Lq^'^ = 3.2 x 10 40 erg/s, which is about
a factor of 3 smaller than the observed luminosity in E > 100 MeV 7-rays L 1 « 10 41 erg/s [90]. To
further constrain the parameters of the model, we evaluate the energy-weighted approximately isotropic
proton flux at 70 EeV. If the source actively emitted UHECRs for at least 70 Myr, from Eq. (73) we
obtain

{E J{E)) = ^)HD{E) dEdt I{x)

w 1.6 x 10 57 eV 3 km" 2 yr" 1 sr" 1 , (77)

in agreement with observations [29]. If we assume circular pixel sizes with 3° radii, the neutrons will be
collected in a pixel representing a solid angle Afi ^ 8.6 x 10" 3 sr. The event rate of (diffuse) protons
coming from the direction of Cen A is found to be

- Q AO /^ 2 /Z?4 7 \ ^ - O HQ nr^ + o/.rr. (J

P 2 dE
S An / (E 4 J) —- = 0.08 events/yr .

d t j El

All in all, in the next 9 yr of operation we expect about 6 direct neutron events against an almost negligible
background. Note that our model also predicts no excess in the direction of M87, in agreement with
observations (see Fig. 8).

We turn now to the discussion of anisotropy. The number of particles with velocity c hitting a unit
area in a unit time in a uniform gas of density n(E, r , t) is n(E, r , t) c. Due to the gradient in the number
density with radial distance from the source, the downward flux at Earth per steradian as a function of
the angle 9 to the source is [133]

J(0, r, t) = n(i Y' t)c (l + a d cos 6) , (79)

where

a d cos6=— — . (80)

n(E,r,t) c

33

The asymmetry parameter can be found by computing the incoming current flux density j = D^V^n
as viewed by an observer on Earth, where r = R + r' . We obtain

n ^0 1 f* dt' -(B?+2R.r'+r' 2 )/4D(t-t') ~( 2R i + 2x 'i)

dEdt (47rD) 3 / 2 J ton (t - t'f/ 2 4D(i - f)

(Ri + x't) r dt'

f

J ton

e -rV[4D(t-t')] (gl)

2(4ttD)3/2 J ton (t-t'f' 2
Near = 0, using the change of variables in Eq. (74), we obtain

tj 2 d£dt (47TD) 3 / 2 4D ^ r 2 y A/x

^ ^° 1 /'(,), (82)

2vr dEdt r 3

where

If 00

l'(x) = -= / du^e~ u .
V 71 " Jl/x

(83)

l/a

Finally, taking R x = R x = and R z = r cos # we obtain [134]

2g(g) /'(*)

For T on oo, < a < 1. Taking our fiducial values E = 70 EeV, 5 nG = 50, and T on = 70, we find
ad = 0.29. This is within ~ la of the anisotropy amplitude ad = d±/ cos«5o = 0.25 ± 0.18 obtained
from the 69 arrival directions, assuming a dipole function for a source model with a maximum value at
CenA: (a ,6 ) = (201.4°, -43.0°).

One caveat is that we assumed that neutrons completely dominate the ultrahigh energy Cen A
emission spectrum; that is

^ = (JVJ + N$)E- S , with JVJ/JVJ«1. (85)

This reduces the number of free parameters in the model. The actual proton-to-neuton fraction depends
on the properties of the source, dominantly by the ratio of photon-to-magnetic energy density. The rela-
tion (85) results from a ~ 0.4, which implies that photopion production - and not proton leakege from
the accelerator region - determines the shape of the cutoff in the spectrum at the source.

EXERCISE 2.2 The assymetry parameter given in Eq. (84) accurately describes the diffuse prop-
agation of CRs if Nq/Nq <C 1. Show that asymmetry parameter from the neutron decay shell is given
by

2D(E) I"(x)

OLd

cr I{x)

where

(1+k) 2 u

(86)

(87)

and k = \{E)/r. Show that in spite of the complicated nature of Eq.(87), the results for ay are very
similar to the ones for the primary diffusion front given above.

34

1

I* %

* *

5 10 15 20

a(hr)

Fig. 15: The nominal arrival directions (a = right ascension, 5 = declination) of SUGAR events with energies
above 4 x 10 10 GeV [138]. Also shown in solid lines are contour maps indicating the circular areas of the celestial
sphere centered at Cen A (indicated by *) with 10° and 25° radii. The dashed lines surrounding several of the
events indicate the uncertainty of each arrival direction; this was found to be about 3° sec 0, where is the shower
zenith angle. From Ref. [81].

An observation demanding of some attention is the distribution of arrival directions of the events
collected by the SUGAR experiment, shown in Fig. 15. We take the data at face value keeping in
mind that the techniques available at the time the experiment was conducted were not as refined as the
techniques currently at our disposal. As one can see in Fig. 15, while the Auger and SUGAR distrbutions
are not dissimilar, the statistics are rather limited and do not support a firm conclusion.

As an attentive reader should know by now, the observed excess in the direction of Cen A can
also be explained using proton directional signals and diffusion of heavy nuclei in a £>-field of about
1 nG [139]. However, if this were the case, a larger anisotropy (as yet unobserved) should be present
at EeV energies [140]. Regarding the preceeding discussion, one should note that the 3° window does
not have an underlying theoretical motivation. Recall that this angular range resulted from a scan of
parameters to maximize signal significance. Cen A covers an eliptical region spanning about 10° along
the major axis; see Fig. 12. Therefore, some care is required to select the region of the sky which is most
likely to maximize the signal-to-noise.

In summary, existing data is consistent with the hypothesis that Cen A dominates the CR sky at
the high end of the spectrum [133]. Auger is in a gifted position to explore Cen A and would provide in
the next 9 yr of operation sufficient statistics to test this hypothesis. The potential detection of neutrons
at Auger can subsequently be validated by the larger aperture of the JEM-EUSO mission.

2.5 Ultrahigh energy cosmic neutrinos
2.5.1 Waxman-Bahcall bound

It is helpful to envision the CR engines as machines where protons are accelerated and (possibly) per-
manently confined by the magnetic fields of the acceleration region. The production of neutrons and
pions and subsequent decay produces neutrinos, 7-rays, and CRs. If the neutrino-emitting source also

35

produces high energy CRs, then pion production must be the principal agent for the high energy cutoff
on the proton spectrum. Conversely, since the protons must undergo sufficient acceleration, inelastic
pion production needs to be small below the cutoff energy; consequently, the plasma must be optically
thin. Since the interaction time for protons is greatly increased over that of neutrons due to magnetic
confinement, the neutrons escape before interacting, and on decay give rise to the observed CR flux. The
foregoing can be summarized as three conditions on the characteristic nucleon interaction time scale 7i nt ;
the neutron decay lifetime r n \ the characteristic cycle time of confinement r cyc i e ; and the total proton
confinement time r conf : (i) r int > r cyc i e ; (ii) r n > r cyc i e ; (Hi) r int < r conf . The first condition
ensures that the protons attain sufficient energy. Conditions (i) and (ii) allow the neutrons to escape the
source before decaying. Condition (Hi) permits sufficient interaction to produce neutrons and neutrinos.
These three conditions together define an optically thin source [141]. A desirable property to reproduce
the almost structureless energy spectrum is that a single type of source will produce cosmic rays with a
smooth spectrum across a wide range of energy.

The cosmic ray flux above the ankle is often summarized as "one 3 x 10 10 GeV particle per
kilometer square per year per steradian." This can be translated into an energy flux [142]

E(EJ CR ) - 3Xl0l ° GeV

(10 10 cm 2 )(3 x 10 7 s)sr
= 10" 7 GeVcm- 2 s- 1 sr- 1 . (88)

From this we can derive the energy density ecR in UHECRs using flux = velocity x density, or

4tt j dE{EJ CR } = ce C R. (89)

This leads to

r-£ max 10 -7 QeV „_ 1Q TeV

4lT f m ' it; .. rv

e C R = — / ~^dE ~ 10" 19 ^3 , (90)

c . h cm z s cm- 3

taking the extreme energies of the accelerator(s) to be E m [ n — 10 10 GeV and £ ma x = 10 12 GeV.
The power required for a population of sources to generate this energy density over the Hubble time

(T H « 10 10 yr) is: e§R°' 10l2] - 5 x 10 44 TeV Mpc~ 3 yr" 1 - 3 x 10 37 erg Mpc~ 3 s" 1 . This works out
to roughly (i) L « 3 x 10 39 erg s _1 per galaxy, (ii) L « 3 x 10 42 erg s _1 per cluster of galaxies, (Hi)
L « 2 x 10 44 erg s _1 per active galaxy, or (iv) J Ldt « 10 52 erg per cosmological GRB [142]. The
coincidence between these numbers and the observed output in electromagnetic energy of these sources
explains why they have emerged as the leading candidates for the CR accelerators.

The energy production rate of protons derived professionally, assuming a cosmological distribu-
tion of sources (with injection spectrum typical of shock acceleration dNo/dE oc E~ 2 ) is [143]

^CR ' 1012 ' ~ 5 x 1044 er g M PC" 3 yr" 1 . (91)

This is within a factor of our back-of-the-envelope estimate (1 TeV =1.6 erg). The energy-dependent
generation rate of CRs is therefore given by

E 2 d ™

.[10 10 ,10 12 ]

dE ln(10 12 /10 10 )

w 10 44 ergMpc" 3 yr- 1 . (92)

The energy density of neutrinos produced through interactions of these protons can be directly
tied to the injection rate of CRs

o dn v 3 ^ o dh

36

where is the fraction of the energy which is injected in protons lost into photopion interactions. The
factor of 3/8 comes from the fact that, close to threshold, roughly half the pions produced are neutral,
thus not generating neutrinos, and one quarter of the energy of charged pion decays goes to electrons
rather than neutrinos. Namely, resonant p~f interactions produce twice as many neutral pions as charged
pions. Direct pion production via virtual meson exchange contributes only about 20% to the total cross
section, but is almost exclusively into 7r + . Hence, p~f interactions produce roughly equal numbers of 7r +
and 7T°. The average neutrino energy from the direct pion decay is {E y ^Y = (1 — r) E n /2 ~ 0.22 E n
and that of the muon is (E^) 71 = (1 + r) ~ 0.78 E^, where r is the ratio of muon to the pion

mass squared. In muon decay, since the has about 1/3 of its parent energy, the average muon neutrino
energy is (E^Y = (1 + r)E n /6 = 0.26 E n .

The "Waxman-Bahcall bound" is defined by the condition = 1

« 2.3 x lO^e^&GeVcm^s^sr" 1 , (94)

where the parameter £ z accounts for the effects of source evolution with redshift, and is expected to
be ^ 3 [144]. For interactions with the ambient gas (i.e. pp rather than p^ collisions), the average
fraction of the total pion energy carried by charged pions is about 2/3, compared to 1/2 in the photopion
channel. In this case, the upper bound given in Eq. (94) is enhanced by 33%. Electron antineutrinos can
also be produced through neutron /3-decay. This contribution, however, turns out to be negligible (see
Appendix A for details).

The actual value of the neutrino flux depends on what fraction of the proton energy is converted to
charged pions (which then decay to neutrinos), i.e. is the ratio of charged pion energy to the emerging
nucleon energy at the source. For resonant photoproduction, the inelasticity is kinematically determined
by requiring equal boosts for the decay products of the A + , giving = E^+ /E n « 0.28, where
and E n are the emerging charged pion and neutron energies, respectively. For pp —> NN + pions, where
N indicates a final state nucleon, the inelasticity is « 0.6 [145]. This then implies that the energy carried
away by charged pions is about equal to the emerging nucleon energy, yielding (with our definition)
« 1.

At production, if all muons decay, the neutrino flux consists of equal fractions of v e , and v^.
Originally, the Waxman-Bahcall bound was presented for the sum of and (neglecting u e ), moti-
vated by the fact that only muon neutrinos are detectable as track events in neutrino telescopes. Since
oscillations in the neutrino sector mix the different species, we chose instead to discuss the sum of all
neutrino flavors. When the effects of oscillations are accounted for, nearly equal numbers of the three
neutrino flavors are expected at Earth [146].

EXERCISE 2.3 The assumption that GRBs are the sources of the observed UHECRs generates
a calculable flux of neutrinos produced when the protons interact with the fireball photons [147]. In
the observer's frame, the spectral photon density (GeV -1 cm -3 ) can be adequately parametrized by a
broken power-law spectrum n^ RB (e 7 ) oc e^, where (3 ~ 1, 2 respectively at energies below and above

e break _ j Me y [148] Show ^

/ inbreak \ — 1

«b(^ > ^ break ) - ID" 13 [Y^GeV ) Cm " 2 S_1 Sr_1 ' (95)

where E^ ieak ~5x 10 5 r| 5 (e^ reak /MeV)- 1 GeV. Recall that in our convention e 7 = Te' T where
is the photon energy measured in the fireball frame. Convince yourself that the non-observation of
extraterrestrial neutrinos from sources other than the Sun and SN 1987a puts the GRB model of UHECR
acceleration on probation [149].

37

If the injected cosmic rays include nuclei heavier than protons, then the neutrino flux expected
from the cosmic ray sources may be modified. Nuclei undergoing acceleration can produce pions, just
as protons do, through interactions with the ambient gas, so the Waxman-Bahcall argument would be
unchanged in this case. However, if interactions with radiation fields dominate over interactions with
matter, the neutrino flux would be suppressed if the cosmic rays are heavy nuclei. This is because the
photodisintegration of nuclei dominates over pion production at all but the very highest energies. Defin-
ing n as the fraction of nuclei heavier than protons in the observed cosmic ray spectrum, the resulting
neutrino flux is then given by E% $^ a11 « (1 — k) E% [150]. The most up-to-date calculation of k

combines a double-fit analysis of the energy and elongation rate measurements to constrain the spectrum
and chemical composition of UHECRs at their sources [151]. Injection models with a wide range of
chemical composition are found to be consistent with observations. In particular, the data are consistent
with a proton dominated spectrum with only a small (1 - 10%) admixture of heavy nuclei. 4

By duplicating the Waxman-Bahcall calculation for Cen A, we obtain an upper limit on the inten-
sity of neutrinos from the direction of the nearest active galaxy,

« 5.0 x 10" 9 GeV cm" 2 s" 1 , (96)

with E v < 10 9 - 5 GeV. For the model introduced in Sec. 2.4, we have 671- ~ 0.28 and therefore a
prediction for the all-flavor neutrino flux

El (/)" c f nA « 1.5 x 10~ 9 GeV cm- 2 s" 1 , (97)

in agreement with the results of [152]. In addition to neutrinos, one also expects a similar flux of photons
at the source, which also carries unperturbed directional information. However, photons at these energies
are difficult to dig out from the hughe proton background.

Although there are no other nearby FRI of this magnitude which can potentially be detected as
point sources, one can integrate over the estimated FRI population out to the horizon to obtain an upper
limit for the diffuse FRI neutrino flux. This quantity is given by [153]

1_ _ . 3

47T

1.5 x 10~ 8 GeV cm- 2 s" 1 sr" 1 , (98)

where TZh — 1 horizon ~ 3 Gpc and n FRI ^8x 10 4 Gpc -3 is the number density [154]. Note that
this flux is about a factor of 3 smaller than the Waxman-Bahcall upper limit. Hence, the reduction in
luminosity of the ensemble of neutrino sources roughly compensates for the presence of distant optically
thin sources whose CR components are hidden by extragalactic magnetic fields.

The diffuse neutrino flux has an additional component originating in the energy losses of UHECRs
en route to Earth [155]. The accumulation of these neutrinos over cosmological time is known as the cos-
mogenic neutrino flux. The GZK reaction chain generating cosmogenic neutrinos is well known [156].
The intermediate state of the reaction P7cmb —> n7T + /pn is dominated by the A + resonance, because
the neutron decay length is smaller than the nucleon mean free path on the CMB. Gamma-rays, pro-
duced via 7T° decay, subsequently cascade electromagnetically on intergalactic radiation fields through
e + e~ pair production followed by inverse Compton scattering. The net result is a pile up of 7-rays at
GeV-TeV energies, just below the threshold for further pair production on the diffuse optical background.
Meanwhile each tt + decays to 3 neutrinos and a positron; the e + readily loses its energy through inverse
Compton scattering on the diffuse radio background or through synchrotron radiation in intergalactic
magnetic fields. As we have seen, the neutrinos carry away about 3/4 of the 7r + energy, therefore the

4 It is important to stress that the essential results of the analysis in Ref. [151] are not altered by the new Auger data [40].

38

energy in cosmogenic neutrinos is about 3/4 of that produced in 7-rays. The intensity of 7-ray pile up
currently provides the most stringent bound on the flux of cosmogenic neutrinos. It is this that we now
turn to study.

2.5.2 Boltzmann equation, universal cosmic ray spectrum, and cosmogenic neutrinos

For a spatially homogeneous distribution of cosmic sources, emitting UHE particles of type i, the comov-
ing number density, Y{(z, E) = rti(z, E)/(l + z) 3 , is governed by a set of one-dimensional (Boltzman)
continuity equations,

Yi = d E (HEYi) + d E {biYi) - V t Y % + ^ J dE 3 ljt Y 3 + Q % ,

(99)

together with the Friedman-Lemaitre equations describing the cosmic expansion rate H(z) as a function
of the redshift z. 5 The first term on the r.h.s. describes adiabatic energy losses [6]. The second term
describes interactions on the cosmic photon backgrounds which can be approximated by continuous
energy losses. The third and fourth terms describe more general interactions involving particle losses
(i -» anything) with interaction rate 1^, and particle generation of the form j -» i. The last term on the
r.h.s., Qi, corresponds to the luminosity density per comoving volume of sources emitting CRs of type i.
We now discuss the calculation of these terms and their scaling with redshift.

The angular-averaged (differential) interaction rate, Ti (7^), appearing on the r.h.s. of Eq. (99) is
defined as

1

Ti(z,Ei) = ^ J d cos 9 J duj{l- (3cos9)n 1 (z,uj)c tot
-1

, (100)

ll3 (z,E u E 3 ) = Ti&Ed^iEitEj), (101)

where n 7 (z,cj) is the energy distribution of background photons at redshift z and dNij/dEj is the
angular- averaged distribution of particles j after interaction. The factor (1 — /3 cos 9) takes into account
the relativistic Doppler shift of the photon density.

In general, any transition i — >► i which can be approximated as 7^ (E, E') ~ 5(E—E f — AE)Ti(E)
with AE/E <C 1 can be replaced in the Boltzmann equations (99) as

-T(E)Y Z (E) + J dE' ju(E f , E) Yi(E f ) -+ d E (b z Y^ , (102)

with bi = AETi ~ — E. The production of electron-positron pairs in the photon background with a
small energy loss is usually approximated as a continuos energy loss (CEL) process. As we have seen
in Sec. 2.3.1, it is also possible to approximate the energy loss in the hadronic cascade due to photopion
production as a CEL

^§(z,E) = bn(z,E) ~ ET p (z,E) - J dE 1 1 E' lpp (z,E,E') , (103)

with &7r(0, E)/E given by (44) and (46). Diffractive processes at high energies with large final state
multiplicities of neutrons and protons ultimately invalidate the CEL approximation. However, the relative
error below 10 12 GeV is less than 15%, so we will use this approximation for a detailed numerical scan

5 This is given by H 2 (z) — Hq [Q m (l + z) 3 + Q\], normalized to its value today of H ~ 100 h kms -1 Mpc -1 , in the
usual "concordance model" dominated by a cosmological constant with Qa ~ 0.7 and a (cold) matter component, Q m ~
0.3 [15]. The time-dependence of the red-shift can be expressed via dz = —dt (1 + z)H.

39

12

E [GeV] red-shift z

Fig. 16: Left: The interaction and decay rates appearing in the Boltzmann equations for the CMB and CIB at
z = 0. Right: Star formation rate (Eq. (Ill) from Ref. [158]) compared with our approximation (110) for the CIB
photon density scaling with redshift. For comparison, the ex (1 + z) 3 scaling of the CMB photon density is also
shown . From Ref. [157].

in the model space of proton spectra. For neutrons with energy less than 10 11 GeV the decay length is
always smaller than the interaction length on the photon backgrounds. It is convenient to approximate
the production of neutrons as ~ Y^ p + T£ n [157]; we have adopted this in all our calculations. In
the left panel of Fig. (16) we show the quantities 6p air /E, 9^6 pair , T p and Ho for comparison.

The redshift scaling of Eqs. (100) and (101) depend on how the photon backgrounds vary with
redshift. The CMB spectral density (GeV -1 cm -3 ) scales adiabatically as:

n 7 (z, u) = (1 + zf n 7 (0, u/(l + z)) , (104)

following from Yy = dE(HEY 7 ) and Y 1 oc a 3 n 7 , where a is the cosmic scale factor. The scaling
behavior of Eq. (104) translates into the following scaling of the quantities Yi and 7^ [21]:

r t (z,E t ) = (i + z) 3 r,(o,(i + z)£;,), (105)

ll3 {z,E u E 3 ) = (l + z)^ ij (0,(l + z)E i ,(l + z)E j ). (106)

The scaling behaviour of b and its derivative c^fr do again depends on the photon background. For the
CMB contribution, we have

b^z.E,) = (l + z) 2 6,(0,(l + z)£;,), (107)
d E b % {z,E % ) = (l + z ) 3 ^(0,(l + z)^). (108)

We now discuss the red-shift scaling of the quantities bi, dsbi, jij and Ti for the case of the cosmic
infrared background. The CIB spectrum has been studied and tabulated for redshifts up to z = 2 [159]
(which we extrapolate slightly to UV energies as seen in Fig. 17). This is consistent with the constraints
on the 7-ray opacity of the universe set by HESS [161], MAGIC [162] and Fermi-LAT [163] from non-
observation of the expected cutoffs in the 7-ray spectra of extragalactic sources. The redshift dependence
is given by:

roc ^

ncmO*, (1 + z)E) = (1 + zf \ dz' ——Q C m(z' , (1 + z')E) , (109)

40

10 4
10 3
10 2

~ 10

I

s 1

, o, 1
0.1

1 io- 2

10" 3

10" 4

10" 5

10" 6
10"

z =

10"

IR/optical (Franceschini et al. '08)
max. radio (Protheroe & Biermann'96)

CMB

10 2

Fig. 17: The energy spectrum of the CMB and the CIB in the IR/optial and radio range at z = 0. The thin dashed
line shows our extrapolation to UV energies. From Ref. [160].

where Qcib is the comoving luminosity density of the sources and we neglect absorption effects other
than expansion. We assume that this follows the star formation rate: Qcm{z,E) oc HsfrX^) Qcib (0, £7).
We can then infer the bolometric evolution to be:

Ncmjz) _ 3 f™dz>Hs F K/(H(z>)(l + z>))

iVciB(O) U + f °°dz>HsFR/(H(z>)(l + z>))' U1U;

where A^cib(^) is the number of infrared-optical photons per proper volume at redshift z. Following a
recent compilation [158, 164], we adopt

(l + zy A Z<1,

7Vi(l + z)-°- 3 l<z<4, (HI)
< 7Vi7V 4 (l + ^)" 3 - 5 z>4.

with appropriate normalization factors, N\ = 2 3 7 and N4 = 5 3 2 (see the right panel of Fig. 16). For
comparison, the CMB evolves as A^cmb(^) /^Vcmb(0) = (1 + z) 3 . To simplify the numerical evaluation
we approximate the evolution with redshift as

n clB (z,u) ~ ^- ^mM nc iB(0,a;/(l + *)) , (112)

and this is shown in the right panel of Fig. 16. The redshift scaling of the quantities 7^ , b{ and

for the CIB is then obtained from the corresponding scaling given for the CMB in Eqs. (105/106) and

(107/108), by multiplying the r.h.s. by a factor N Qm {z) /Ncib(0)/(1 + zf .

We have little direct knowledge of the cosmic radio background. An estimate made using the RAE
satellite [165] is often used to calculate the cascading of UHE photons. A theoretical estimate has been
made [166] of the intensity down to kHz frequencies, based on the observed luminosity function and
radio spectra of normal galaxies and radio galaxies, although there are large uncertainties in the assumed
evolution. The calculated values are about a factor of ^ 2 above the measurements and to ensure maximal

41

energy transfer in the cascade we will adopt this estimate and assume the same redshift scaling as the
cosmic infrared/optical background. We summarize the adopted cosmic radiation backgrounds in Fig. 17.

The emission rate of CR protons per comoving volume is assumed, as per usual practice, to follow
a power-law:

{/_ [E j i?min) E < E m [ n ,
1 £min<£<£max, (H3)

f+(E/ £ max ) ^max < E .

We will consider spectral indices 7 in the range 2-1-3. The functions f±(x) = x ±2 exp(l — x ±2 ) in
Eq. (113) smoothly turn off the contribution below i? m in and above i? m ax- We set £ max = 10 12 GeV
in the following and vary E m i n in the range 10 8 5 -7- 10 10 GeV, corresponding to a galactic-extragalactic
crossover between the second knee and the ankle in the CR spectrum.

The cosmic evolution of the spectral emission rate per comoving volume is parameterized as:

Qp(z,E)=H(z)Qp(0,E). (114)
For simplicity, we use the standard approximation

n(z) = (l + z) n &(z m ^-z), (115)

with z max = 2. Note that the dilution of the source density due to the Hubble expansion is taken care of
since Q is the comoving density, i.e. for no evolution we would simply have H = 1. We consider cosmic
evolution of UHECR sources with n in the range 2-1-6.

We can express the system of partial integro-differential equations (99) in terms of a system of
ordinary integro-differential equations [157],

£ i = -HS i -b i (z,S i ), (116)
Zi = [Pi(z,£i) - Ti(z,£i)] Z i + (l + z)Qf{z,£ i ) , (117)

where we have defined f3i(z,E) = dEbi(z,E) and Zi(z,E) = (1 + z)Yi(z,£i(z,E)). The quantity
£i(z, E) gives the energy that a particle of type i had at redshift z if we observe it today with energy E
and take into account CEL. The effective source term in Eq. (1 17) is

Qf{z,£ i {z,E)) = Q i + Y, J dEd E £ 3l3% {z,£ 3 ,E % )^-, (118)

3

where £j (z , E) and Zj (z, E) are subject to the boundary conditions £j (0, E) = E and Zj (^ m ax, E) = 0.
The flux at z = can be expressed as

ME) = ^Zi(0,E)

^ />oo r pz

— — dzexr> / dz'
4vr J [Jo

j d E b i (z',e i (z',E))-T i (z',ei(z',E))
(l + z')H(z>)

H(z)

For the numerical evaluation of the Boltzmann equations (99) it is convenient to solve [160]

-H(z)(l + z)d z 3fi(z, E) = -T(z, £)m(z, E) + -^—^d E [b(z, £)%{z, E)\ + (1 + z)Qf(z, S) .

(120)

where we have defined £ — (l+z)E, and 3?i(z, E) = (l+z)Yi(z, $), subject to the boundary condition
^K^max, E) = 0. The effective source term becomes Qf ff (z, <§) = Q% + J2j J dEj jji(z, $j,Ei) 3?j.
In our calculation we use logarithmic energy bins with size Alog 10 E = 0.05 between 10 5 GeV and

42

10 5 GeV. For the numerical evolution of the differential Eq. (1 16) we use a simple Euler method with a
step-size Az = 1CT 4 . The corresponding step-size in the propagation distance Ar = cAt is then always
smaller than the proton interaction length. The flux at z = is given by Ji{E) = i^(0, E)/ (47r).

Electromagnetic (EM) interactions of photons and leptons with the extra-galactic background light
and magnetic field can happen on time-scales much shorter than their production rates. It is convenient
to account for these contributions during the proton propagation as fast developing electro-magnetic cas-
cades at a fixed redshift. We will use the efficient method of "matrix doubling" [167] for the calculation
of the cascades. Since the cascade 7-ray flux is mainly in the GeV-TeV region and has an almost univer-
sal shape here, it is numerically much more efficient to calculate the total energy density cj cas injected
into the cascade and compare this value to the limit imposed by Fermi-LAT. The total energy density (eV
cm -3 ) of EM radiation from proton propagation in the past is given as

^ cas = J dEEncss^E) = JdtJ dE b ™^' z ^ np(z,E), (121)

where n p (z, E) is the physical energy density at redshift z, defined via n p (z, E) = (1 + z) 3 Y p (z, E)\
details on the calculation are given in Appendix B. The continuous energy loss of protons into the
cascade, denoted by 6 cas , is in the form of electron, positron and 7-ray production in BH (&bh) an d
photopion (b n ) interactions. We derive the BH and photopion contribution to cj cas separately. For the
photopion contribution we estimate

b n (z,E) ~ J dE'E' [ lpe -(z,E,E') + lpe+ (z,E,E') + lpi (z,E,E')] , (122)

where the angular-averaged distribution of particle j after the interaction are determined using the Monte
Carlo package SOPHIA [168]. For the energy loss via BH pair production we use (37). Note, that since
the photopion contribution in the cascade is dominated at the GZK cutoff, its contribution should increase
relative to BH pair production with increasing crossover energy and, hence, also the associated neutrino
fluxes after normalization to 7-ray and CR data.

Next, we study the constraint set by the diffuse 7-ray background on all-proton models of ex-
tragalactic CRs. We parametrize our ignorance of the crossover energy - which marks the transition
between the galactic and extragalactic components - as a variable low energy cutoff in the proton in-
jection rate. By fitting only to CR data above the crossover energy, taken to be between 10 8 5 GeV and
10 10 GeV, we determine the statistically preferred values of the spectral index 7, cosmic source density
evolution index n, and crossover energy E m [ n by a goodness-of-fit (GOF) test of the HiRes data, taking
into account the energy resolution of about 25%. For each model we check that the total energy density
of the EM cascade is below a critical value inferred from the recent measurement of the extragalactic
7-ray background by the Fermi-LAT Collaboration [169].

Given the acceptance A{ (in units of area per unit time per unit solid angle) of the experiment for
the energy bin i centered at E{ and with bin width A^, and the energy scale uncertainty of the experiment,
(Te s the number of expected events in the bin is given by

Ni{n^N^S) = Ai / J p Nqu (E)dE, (123)

where

J P N ,n,,(E) = ^n p (0,E) (124)

is the proton flux arriving at the detector corresponding to a proton source luminosity as in Eq. (113), with
the cosmic evolution of the source density given by Eqs. (114) and (115). The parameter S in Eq. (123)
above is a fractional energy-scale shift that reflects the energy-scale uncertainty of the experiment, and
Nq is the normalization of the proton source luminosity.

43

n n n n

Fig. 18: Left Panel: Goodness-of-fit test of the HiRes data [28]. We show the 68% (pink), 95% (blue) and 99%
(magenta) confidence levels of the injection index 7 and the cosmic evolution index n . The black lines indicate the
allowed regions before the cascade (co> cas ) bound is imposed. Right Panel: The corresponding energy density in
the EM cascade. From Ref. [160].

The probability distribution of events in the i-th bin is of the Poisson form with mean A^. Corre-
spondingly the r-dimensional (r being the number of bins of the experiment with Ei > E^^) probability
distribution for a set of non-negative integer numbers k = {fci, ...k r }, Pj^{n, 7, iVo, 5), is just the product
of the individual Poisson distributions. According to this r-dimensional probability distribution, the ex-
perimental result iV exp = {A^ xp , A^ xp } has a probability Pfi exp (n, 7, iVo, S) and, correspondingly,
the experimental probability after marginalizing over the energy scale uncertainty and normalization is:

(n,7,JV ,<5), (125)

where the maximization is made within some prior for 6 and Nq. For the energy shift 6 we have used
two forms for the prior, either a top-hat spanning the energy- scale uncertainty of the experiment, cfe s , or
a gaussian prior of width <je s •

For No we impose the prior arising from requiring consistency with the Fermi-LAT measure-
ments [169] of the diffuse extra-galactic 7-ray background. In order to do so, we obtain the total energy
density of EM radiation from the proton propagation using Eq. (121) and we require following Ref. [170]:

,(N 0l n, 7) < 5.8 x 10" 7 eV/cm 3 . (126)

The marginalization in Eq. (125) also determines iVQ est and S hest for the model, which are the
values of the energy shift and normalization that yield the best description of the experimental CR data,
subject to the constraint imposed by the Fermi-LAT measurement.

Altogether the model is compatible with the experimental results at given GOF if

J2 Pk(n, 7, < eSt , S best )0 \p n (n, 7, iV bes \ 5 best ) - P exp (n, 7)

< GOF . (127)

Technically, this is computed by generating a large number iV re p of replica experiments according to
the probability distribution P%(n, 7, NQ est , S hest ) and imposing the fraction F of those which satisfy

P^n, 7, JV , 6 hest ) > P exp (n, 7) to be F < GOF.

44

GOF99%

>

E min =10 17 ' 5 eV

* HiRes II

• HiRes I
o Auger

G0F99%

>
o

E mIn =10 18 eV

* HiRes II
■ HiRes I
o Auger

E(GeV)

E(GeV)

Fig. 19: The allowed proton flux (at the 99% confidence level) for increasing crossover energy E m { n . Each fit of
the proton spectrum is marginalized with respect to the experimental energy uncertainty and we show the shifted
predictions in comparison to the HiRes central values [28]. For comparison we also show the Auger data [29, 30]
which has not been included in the fit. From Ref. [160].

With this method we determine the value of (n, 7) parameters that are compatible with the HiRes
I and HiRes II experiments [28]. In the left panel of Fig. 18, we plot the regions with GOF 64%,
95% and 99% for four values of the minimum (i.e. crossover) energy. In the right panel, we show
the corresponding ranges of cj cas? best f° r the models as a function of the cosmic evolution index n. In
order to display explicitly the impact of the constraint from the Fermi-LAT measurements of the diffuse
extra-galactic 7-ray background (126), we show the corresponding GOF regions without imposing that
constraint. In Table. 3 we list the parameters corresponding to the best-fit models and to the models
with minimal and maximal contributions to uo^ and cj cas = uo^ + c^bh at the 99% C.L., together with
the corresponding energy shifts which give best fits to the HiRes I and HiRes II data. We also show the
parameters for the models with maximum uj^ and cj cas without imposition of the Fermi-LAT constraint.

It is interesting to compare the allowed values of n with those in the cosmological evolution of
UHECR candidate sources. For GRBs [171],

^grb(^) = (1 + ^) L4 ^sfrW, (128)

for AGNs [172,173]

'HagnO)

(l + z) 5 - u z < 1.7,

N lm7 1.7<z<2.7, (129)

45

n

Fig. 20: Systematic effect of the experimental energy resolution on the fitted spectral index 7 and cosmological
evolution parameter n. For illustration we show the dependence of the 95% C.L. bound for a crossover energy of
10 18 eV. The blue contour corresponds to the region shown in Fig. 18 assuming an uncorrelated energy shift of
25% in both data sets (HiRes I and II) [28], for a flat prior ("top-hat" distribution). The red dashed curve assumes
correlated errors of the energy resolution in both data sets. The black dotted curve shows the result for uncorrelated
errors with a Gaussian prior, and the dashed-dotted line shows uncorrelated errors with a flat prior, but with a lower
uncertainty of 15%. From Ref. [160].

with TVi.r 2.7 5 and N 2 . 7

10 043 , and for QSO [174]

f(l + z) 3
= < (1 + 1.9) 3

[(l + 1.9) 3 exp{(2.7-z)/2.7},

z < 1.9,

1.9 < z < 2.7, (130)
z > 2.7.

As an illustration of the agreement with the CR data, in Fig. 18 we show the range of proton
fluxes corresponding to models with GOF 99% or better for increasing crossover energies E m [ n . As
discussed above, each fit of the proton spectra is marginalized with respect to the experimental energy
scale uncertainty and we show the shifted predictions with S hest in comparison to the HiRes data at
central value. In the figure we also show the results from Auger [29, 30], though these have not been
included in the analysis.

These results are obtained assuming an energy scale uncertainty <je s = 5% with a "top-hat" prior
for the corresponding energy shifts which are taken to be uncorrelated for HiRes I and HiRes II. In Fig. 20
we explore the dependence of the results on these assumptions by using a different form for the prior,
assuming the energy shifts to be correlated between the two experiments, or reducing the uncertainty to
ce s = 15%. As seen in the figure, the main effect is associated with the reduction of the energy scale
uncertainty which, as expected, results in a worsening of the GOF for models with larger n. This is
directly related to the normalization constraint from Eq. (126). If one naively ignores the energy scale
uncertainty, the constraint in Eq. (126) rules out models with n > 3 (the precise value depending on the
assumed E m { n ). However, once the energy scale uncertainty is included, the constraint of Eq. (126) plays
a weaker role on the determination of the GOF of the models. It does, however, imply a maximum value
of A^Q est which, as we will see, impacts the corresponding ranges of neutrino fluxes.

46

Table 3: Cosmic ray source parameters which best fit the HiRes data [28], along with those which yield minimal
and maximal contributions to uj^ (i.e. neutrino fluxes) and uj cas = uj^ -f co>bh (i.e. 7-ray fluxes), all at the 99% C.L.

£ min = 10 8 - 5 GeV

E m[n = 10 9 GeV

mouei

n

7

1 > a

^cas

^/best ^J/best

n 7

, , a
^cas

^/best

^7/best

fit with Fermi-LAT bound:

best fit

3.50

2.49

5.8

0.005 0.

3.20 2.52

5.2

0.050

0.045

min. uj cas

4.50

2.31

4.4

-0.235 -0.245

2.25 2.47

1.7

-0.120

-0.150

max. u cas

4.60

2.36

5.8

-0.185 -0.175

3.35 2.55

5.8

0.050

0.060

min. uj^

2.00

2.67

4.9

0.215 0.235

2.00 2.51

1.8

-0.070

-0.095

max. uj n

4.80

2.29

5.8

-0.220 -0.215

5.10 2.29

5.8

-0.250

-0.250

fit without Fermi-LAT bound:

max. uj cas

4.45

2.44

15

0.135 0.155

5.25 2.36

27

0.205

0.205

max.

4.80

2.36

14

0.050 0.055

5.30 2.35

26

0.190

0.190

E min = 10 9 - 5 GeV

E min = 10 10 GeV

model

n

7

, , a
^cas

^/best ^//best

n 7

1 > a

^cas

^/best

dllbest

fit with Fermi-LAT bound:

best fit

4.05

2.47

5.8

0.015 0.005

4.60 2.50

4.4

-0.030

-0.065

min. uj cas

2.00

2.45

1.4

-0.050 -0.060

2.00 2.88

0.44

-0.220

-0.250

max. u cas

4.95

2.37

5.8

-0.165 -0.160

4.45 2.13

5.8

0.130

0.090

min.

2.00

2.63

2.1

0.075 0.070

2.00 2.88

0.44

-0.220

-0.250

max. uj n

5.35

2.28

5.8

-0.240 -0.250

4.40 2.10

5.8

0.145

0.100

fit without Fermi-LAT bound:

max. uj cas

6.00

2.49

30

0.120 0.135

6.00 2.14

23

0.250

0.210

max.

6.00

2.47

29

0.120 0.125

6.00 2.10

23

0.250

0.210

a in units of 10" 7 eV/cm 3

The corresponding range of 7-ray and cosmogenic neutrino fluxes (summed over flavor) is shown
in Fig. 21 for models with minimal and maximal energy density at the 99% C.L. As expected, the
maximum 7-ray fluxes are consistent with the Fermi-LAT data within the errors. For illustration, we also
show as a dotted line the "naive" 7-ray limit E 2 J C8iS < ccj™| x /47rlog(TeV/GeV), corresponding to a
7-ray flux in the GeV-TeV range which saturates the energy density (126).

The cosmogenic neutrino and photon fluxes are also sensitive to the primary composition. For
example, if the primaries are iron nuclei one would expect a considerably lower flux on both photons and
neutrinos, rendering photon pile-up measurements less helpful in constraining cosmogenic fluxes. (See
Appendix A for details).

2.5.3 Upper limits on the cosmic neutrino flux

High energy neutrino detection is one of the experimental challenges in particle astrophysics for the
forthcoming years. It is widely believed that one of the most appropriate techniques for neutrino detection
consists of detecting the Cherenkov light from muons or showers produced by the neutrino interactions
in underground water or ice. For a recent review see e.g. [175]. This allows instrumentation of large
enough volumes to compensate for both the low neutrino cross section and the low fluxes expected.

47

i — i HiRes I&II

i — i Fermi LAT

— p (best fit)

— u,u (99% C.L.)

— 7 (99% C.L.)

^min = 10 19 eV

1 1 maximal cascade

I

\

/IS \

- i — i HiRes I&II \\ /

i — i Fermi LAT \\ /I

— P (best fit) \V, //

r — u,u (99% C.L.) Sk / /

— 7 (99% C.L.) \\ / /

0.1 1 10 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 10 10 10 11 10 1 '
E [GeV]

0.1 1 10 10 2 10 3 10 4 10 5 10 6 10 7 10 8

E [GeV]

10 10 10 11 10 12

Fig. 21: Comparison of proton, neutrino and 7-ray fluxes for different crossover energies. We show the best- fit
values (solid lines) as well as neutrino and 7-ray fluxes within the 99% C.L. with minimal and maximal energy
density (dashed lines). The values of the corresponding model parameters can be found in Table. 3. The dotted line
labeled "maximal cascade" indicates the approximate limit E 2 J cas < cu^f^/Air log(TeV/GeV), corresponding
to a 7-ray flux in the GeV-TeV range saturating the energy density (126). The 7-ray fluxes are marginally consis-
tent at the 99% C.L. with the highest energy measurements by Fermi-LAT. The contribution around 100 GeV is
somewhat uncertain due to uncertainties in the cosmic infrared background. From Ref. [160].

There are several projects under way to build sufficiently large detectors to measure the expected signals
from a variety of neutrino sources. The IceCube facility, deployed near the Amundsen-Scott station, is
the largest neutrino telescope in the world [176]. It comprises a cubic-kilometer of ultra-clear ice about
a mile below the South Pole surface, instrumented with long strings of sensitive photon detectors which
record light produced when neutrinos interact in the ice.

CR experiments, like Auger, provide a complementary technique for UHECv detection by search-
ing for deeply-developing, large zenith angle (> 75°) showers [177]. At these large angles, hadron-
induced showers traverse the equivalent of several atmospheres before reaching detectors at the ground.
Beyond about 2 atmospheres, most of the electromagnetic component of a shower is extinguished and
only very high energy muons survive. Consequently, a hadron-induced shower front is relatively flat and
the shower particles arrive within a narrow time window (Fig. 22, top panel). In contrast, a neutrino
shower exhibits characteristics similar to those of a vertical shower, which has a more curved front and a
wider distribution in particle arrival times due to the large number of lower energy electrons and photons.
Furthermore, the "early" part of the shower will tend to be dominated by the electromagnetic component,
while "late" portion will be enriched with tightly bunched muons (Fig. 22, middle panel). Using these
characteristic features, it is possible to distinguish neutrino induced events from background hadronic

48

showers. Moreover, because of full flavor mixing, tau neutrinos are expected to be as abundant as other
species in the cosmic flux. Tau neutrinos can interact in the Earth's crust, producing r leptons which may
decay above the Auger detectors [178-180] (Fig. 22, bottom panel). Details on how such events can be
selected at the Auger Observatory are discussed in [181, 182].

Fig. 22: Schematic illustration of the properties of a hadron-induced shower (top), an z/-induced nearly horizontal
shower (middle) and a ^-induced earth skimming shower (bottom). Note that only up-going showers resulting
from t neutrino interactions in the Earth can be detected with any efficiency using the surface array. In contrast,
all three neutrino species can be detected in down-going showers. Also note from the inset in the lower panel, that
the incident r can experience several CC interactions and decays and thereby undergo a regeneration process. The
Andes mountain range lies to the west of the observatory, and provides roughly an additional 20% target volume
for v T interactions.

So far, no neutrino candidates have been observed resulting in upper limits on the diffuse flux of
neutrinos. For the case of up-going r neutrinos in the energy range 2 x 10 8 GeV < E v < 2 x 10 10 GeV,
assuming a diffuse spectrum of the form E~ 2 , the current 90%CL bound, is [183]

El ~ 4.7;^ x 10~ 8 GeV cm- 2 s" 1 sr" 1 . (131)

Though Auger was not designed specifically as a neutrino detector, it is interesting to note that it exhibits
good sensitivity in an energy regime complementary to those available to other dedicated instruments.
For example, the 90% CL upper limit from IceCube, for neutrinos of all flavors in the energy range

49

2.0 x 10 6 GeV < E v < 6.3 x 10 9 GeV, is [184]

E 2 ^i/aii _ 3 g x 1Q -8 Q e y cm -2 s -l gr -l _ ( j 32)

We now derive model-independent bounds on the total neutrino flux. The event rate for quasi-
horizontal deep showers is

N = E / dEi Na °iN^x(Ei) £{E X ) , (133)

i,x J

where the sum is over all neutrino species i = v e , v e , v^, v^, v T , v T , and all final states X. Na =
6.022 x 10 23 is Avogadro's number, and <fr l is the source flux of neutrino species i, a as usual denotes
the cross section, and £ is the exposure measured in cm 3 w.e. sr time. The Pierre Auger Collaboration
has searched for quasi-horizontal showers that are deeply-penetrating [183]. There are no events that
unambiguously pass all the experimental cuts, with zero events expected from hadronic background.
This implies an upper bound of 2.4 events at 90%CL from neutrino fluxes [185]. Note that if the number
of events integrated over energy is bounded by 2.4, then it is certainly true bin by bin in energy. Thus,
using Eq. (133) one obtains

V / dE i N A <S> i (E i )a iN ^ x (E i )£(E i ) <2.4, (134)

7^ Ja

at 90% CL for some interval A. Here, the sum over X takes into account charge and neutral current
processes. In a logarithmic interval A where a single power law approximation

&(Ei) (JiN^xiEi) £(Ei) ~ Ef (135)

is valid, a straightforward calculation shows that

[ {E)e ' ^E^ a iN ^ x £ = (a iN ^ x £ E % **) ^ A , (136)

J(E)e- A / 2 Ei

where S = (a + l)A/2 and (A) denotes the quantity A evaluated at the center of the logarithmic
interval [186]. The parameter a = 0.363 + f3 — 7, where 0.363 is the power law index of the SM
neutrino cross section [187] and f3 and —7 are the power law indices (in the interval A) of the exposure
and flux <3> z , respectively. Since sinhS/S > 1, a conservative bound may be obtained from Eqs. (134)
and (136):

N A ^2((TiN^x(Ei)) (S(Ei)) (Ei&) < 2.4/A . (137)

i,X

By taking A = 1 as a likely interval in which the single power law behavior is valid (this corresponds
to one e-folding of energy), it is straightforward to obtain upper limits on the neutrino flux. The model-
independent upper limits on the total neutrino flux, derived using an equivalent of 0.8 yr of full Auger
exposure, are collected in Table 4 [183].

3 Phenomenology of extensive air showers

3.1 Systematic uncertainties in air shower measurements from hadronic interaction models

Uncertainties in hadronic interactions at ultrahigh energies constitute one of the most problematic sources
of systematic error in the analysis of air showers. This section will explain the two principal schools of
thought for extrapolating collider data to ultrahigh energies.

Soft multiparticle production with small transverse momenta with respect to the collision axis is
a dominant feature of most hadronic events at cm. energies 10 GeV < < 62 GeV (see e.g., [188,

50

Table 4: Model-independent upper limits on the neutrino flux at 90% CL.

E v (GeV) (Ey gjg (cm" 2 sr" 1 s" 1 )

1 x 10 8 4.3 x 10- 14

3 x 10 8 5.3 x 1(T 15

1 x 10 9 1.2 x l(r 15

3 x 10 9 4.7 x l(r 16

1 x 10 10 2.2 x 10~ 16

3 x 10 10 1.3 x 10- 16

1 x 10 11 7.2 x 1CT 17

3 x 10 11 4.3 x 1(T 17

189]). Despite the fact that strict calculations based on ordinary QCD perturbation theory are not feasible,
there are some phenomenological models that successfully take into account the main properties of the
soft diffractive processes. These models, inspired by 1/N QCD expansion are also supplemented with
generally accepted theoretical principles like duality, unitarity, Regge behavior, and parton structure. The
interactions are no longer described by single particle exchange, but by highly complicated modes known
as Reggeons. Up to about 62 GeV, the slow growth of the cross section with y/s is driven by a dominant
contribution of a special Reggeon, the Pomeron.

At higher energies, semihard (SH) interactions arising from the hard scattering of partons that
carry only a very small fraction of the momenta of their parent hadrons can compete successfully with
soft processes [190-197]. These semihard interactions lead to the "minijet" phenomenon, i.e. jets with
transverse energy (Et = \p T \) much smaller than the total cm. energy. Such low-p T processes cannot
be identified by jet finding algorithms, but (unlike soft processes) still they can be calculated using
perturbative QCD. The cross section for SH interactions is described by

"QCD(«,j£ in ) = £ / — [— d\t\ si/iOri, \t\) x 2 f j (x 2 , \t\) , (138)

where x\ and X2 are the fractions of the momenta of the parent hadrons carried by the partons which col-
lide, d&ij/d\i\ is the cross section for scattering of partons of types i and j according to elementary QCD
diagrams, fi and fj are parton distribution functions (PDFs), s = x\ X28 and —i=s (1— cos #*)/2 = Q 2
are the Mandelstam variables for this parton-parton process, and the sum is over all parton species. Here,

p T = E^ sintf jet = ^ sinr , (139)

and

p„ =£] e a t b COS^jet, (140)

where E 1 ^ is the energy of the jet in the lab frame, $j et the angle of the jet with respect to the beam
direction in the lab frame, and #* is the angle of the jet with respect to the beam direction in the cm.
frame of the elastic parton-parton collision. This implies that for small p 2 T « Q 2 . The integration
limits satisfy

Q 2 min < \t\ <s/2, (141)

where Qmin = 1 — 2 GeV is the minimal momentum transfer. The measured minijet cross sections indi-
cate that the onset of SH interactions has just occurred by CERN SPS energies (y/s > 200 GeV [198].

A first source of uncertainty in modeling cosmic ray interactions at ultrahigh energy is encoded in
the extrapolation of the measured parton densities several orders of magnitude down to low x. Primary

51

protons that impact on the upper atmosphere with energy ~ 10 11 GeV yield partons with x = 2p* / ~
vn^j ~ 10 -7 , whereas current data on quark and gluon densities are only available for x > 10 -4 to
within an experimental accuracy of 3% for Q 2 « 20 GeV 2 [199]. In Fig. 23 we show the region of the
x — Q 2 plane probed by HI, ZEUS, and fixed target experiments. In addition, extrapolation of HERA
data to UHECR interactions assumes universality of the PDFs. This assumption, based on the QCD
factorization conjecture (the cross section of Eq. (138) can always be written in a form which factorizes
the parton densities and the hard interaction processes irrespective of the order in perturbation theory
and the particular hard process) holds in the limit Q 2 » Aqcd, where Aqcd ~ 200 MeV is the QCD
renormalization scale.

For large Q 2 and not too small x, the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equa-
tions [201-204]

9 (g^Q 2 )\ = a,(Q 2 ) fP qq P q9 \ - (q(x,Q 2 )\ (U2i

d\nQ 2 \g(x,Q 2 )J 2n \P gq Pj ™ \g(x,Q 2 )J K }

successfully predict the Q 2 dependence of the quark and gluon densities (q and respectively). Here,
a s = g 2 /(4:7r), with g the strong coupling constant. The splitting functions indicate the probability
of finding a daughter parton i in the parent parton j with a given fraction of parton j momentum. This
probability will depend on the number of splittings allowed in the approximation. In the double-leading-
logarithmic approximation, that is lim^^o ln(l/x) and limg2^ 00 1ii((5 2 /Aqcd), the DGLAP equations
predict a steeply rising gluon density, xg ~ x -cu , which dominates the quark density at low x. This
prediction is in agreement with the experimental results from HERA shown in Fig. 24 [206,207]. HERA
data are found to be consistent with a power law, xg(x, Q 2 ) ~ x~ Ali , with an exponent Ah between 0.3
and 0.4 [208].

52

Q! 20

WD

* 17.5

15
12.5
10
7.5
5
2.5

Q 2 =20 GeV 2

Q 2 =200 GeV 2

HI NLO-QCD Fit 2000

total uncert.
exp. uncert.

ZEUS NLO-QCD Fit
(Prel.) 2001

exp. uncert.

10

Fig. 24: Gluon momentum distributions xg(x, Q 2 ) in the proton as measured by the ZEUS and HI experiments at
various Q 2 . From Ref. [205].

The high energy minijet cross section is then determined by the dominant gluon distribution

^qcd^pT) « f /2 d\t\ ^| x ig ( Xl ,\i\) x 2 g(x 2 ,\t\) , (143)

J xi J x 2 Jq2. d\t\

where the integration limits satisfy

x\x 2 s > 2\t\ > 2Q

2

min '

(144)

Furthermore, because da/d\i\ is peaked at the low end of the \i\ integration (see e.g [209]), the high
energy behavior of ctqcd is controlled (via the lower limits of the x\, x 2 integrations) by the small-x
behavior of the gluons [210]

( \ f dx ^ -A H f dx 2 -A H A h i

J2Q 2 . Is x l J2Q 2 . Ixi8 x< 2

s

(145)

53

This estimate is, of course, too simplistic. At sufficiently small x, gluon shadowing corrections suppress
the singular x~ Ah behavior of xg and hence suppress the power growth of ctqcd with increasing s.

Although we have shown that the onset of semihard processes is an unambiguous prediction of
QCD, in practice it is difficult to isolate these contributions from the soft interactions. Experimental
evidence indicates that SH interactions can essentially be neglected up to and throughout the CERN ISR
energy regime, y/s < 62 GeV. Therefore, measurements made in this energy region can be used to
model the soft interactions. A reasonable approach introduced in [21 1] is to base the extrapolation of the
soft interactions on the assumption of geometrical scaling [212], which is observed to be true throughout
the ISR energy range [213,214]. To this end, we introduce the standard partial- wave amplitude in impact-
parameter space f(s,b), which is the Fourier transform of the elastic pp (or pp) scattering amplitude. (We
neglect any difference between pp and pp for > 200 GeV.) Geometrical scaling (GS) corresponds
to the assumption that /, which a priori is a function of two dimensional variables b and s, depends only
upon one dimensional variable (3 = b/R(s), where R is the energy dependent radius of the proton, i.e.

f(s,b) = f GS (/3 = b/R(s)). (146)

Physically, this means that the opaqueness of the proton remains constant with rising energy and that the
increase of the total cross section, <T to t> in the ISR energy range reflects a steady growth of the radius
R(s). An immediate obvious consequence of GS is that the partial wave at b = should be independent
of energy

/M = 0) = /gs(/? = 0). (147)

Another consequence is that the ratio of elastic scattering to total cross section, a e \(s)/a to t(s), should
be energy independent. This follows from

otot = 8tt J lmf(s,b)bdb

= 8tvR\s) Jlmf GS (P)Pdp, (148)

GS

and

<Tel = 8^ J \f(s,b)\ 2 bdb

= 8ttR 2 (s) J |/ GS (/3)| 2 /3d/3. (149)

GS

To determine the gross features at high energies we can assume that the elastic amplitude has a simple
form

F(s,t) = ia tot (s) e Bt /\ (150)

with B the slope parameter that measures the size of the proton [215]. This is a reasonable assumption:
the amplitude is predominantly imaginary, and the exponential behavior observed for \t\ < 0.5 GeV 2
gives the bulk of the elastic cross section. Now, the Fourier transform f(s,b) of the elastic amplitude
t) given by (150) has a Gaussian form in impact parameter space

/M) = ^, (151)

and it follows that

Im/( S ,6 = 0) = ^| = 2 ^. (152)

87T.O cr tot

54

Equation (152) offers a very clear way to see the breakdown of GS and to identify semihard interactions
from the growth of the central partial wave.

In general unitarity requires Im/(s, 6) < \, which in turm implies cr e \/(Jtot < \ [215]. This seems
to indicate that the Gaussian form (151) may not longer be applicable at ultrahigh energies, but rather it
is expected that the proton will approximate a "black disk" of radius bo, i.e. /(s, b) = | for < b < bo
and zero for b>bo. Then a e \ ~ ^cr to t — ^b\.

In order to satisfy the unitarity constraints, it is convenient to introduce the eikonal x defined by

f( s ,b)= t -{l-exp[i x (s,b)]} , (153)

where Imx > 0. If we neglect for the moment the shadowing corrections to the PDFs and take xg oc
X -A H j n ^ sma n_ x ii m it W e obtain, as explained above, power growth of the cross section for SH
interactions, oqcd ~ s Au and Imx(s,6 = 0) » 1 as s — >► oc. Indeed we expect Imx ^ 1 ( an d
unitarity to be saturated) for a range of b about 6 = 0. Then we have

roc

<rtot = 4^ / bdbQ(b -b)
Jo

rbo(s)

~ 4tt / bdb = 2irbl, (154)
Jo

with x — Xsh an( i where bo(s) is such that

Mo(*))^l. (155)

Hereafter, we ignore the small real part of the scattering amplitude, which is good approximation
at high energies. The unitarized elastic, inelastic, and total cross sections (considering now a real eikonal
function) are given by [216-219]

cr e i = 27r J dbb {1 -exp[-x soft (5,6) -x SH ( 5 ^)]} 2 > ( 156 )
= 2n J dbb{l-exp[-2x soit (s,b)-2x Sii (s,b)}} , (157)
cr t ot = 47r y dbb {1 -exp[-x soft (s,6) -X SH ( 5 ^)]} ?

Cinel

(158)

where the scattering is compounded as a sum of QCD ladders via SH and soft processes through the
eikonals x SH and x soft -

Now, if the eikonal function, x(s, b) = x soft (<§, 6) + Xsh ( s ^) = indicates the mean number
of partonic interaction pairs at impact parameter 6, the probability p n for having n independent partonic
interactions using Poisson statistics reads, p n = (X n /n\) e~ x . Therefore, the factor 1 — e~ 2x = Yl^Li Pn
in Eq. (157) can be interpreted semiclassically as the probability that at least 1 of the 2 protons is broken
up in a collision at impact parameter b. With this in mind, the inelastic cross section is simply the
integral over all collision impact parameters of the probability of having at least 1 interaction, yielding
a mean minijet multiplicity of (n m i n ij et ) ~ CQCD/^inei [220]. The leading contenders to approximate
the (unknown) cross sections at cosmic ray energies, SIBYLL [221] and QGSJET [222], share the eikonal
approximation but differ in their ansdtse for the eikonals. In both cases, the core of dominant scattering
at very high energies is the SH cross section given in Eq. (138),

X S H = l*QCD(s,p™ in )A(s,b), (159)

55

where the normalized profile function, 2tt dbbA{s 1 b) = 1, indicates the distribution of partons in
the plane transverse to the collision axis.

In the QGS JET-like models, the core of the SH eikonal is dressed with a soft-pomeron pre-evolution
factor. This amounts to taking a parton distribution which is Gaussian in the transverse coordinate dis-
tance 6,

e -b 2 /R 2 (s)

A(s,b)= , (160)

7rR Z {S)

with R being a parameter. For a QCD cross section dependence, oqcd ~ s Ah > one gets for a Gaussian
profile

b 2 (s) -i? 2 A H Ins (161)

and at high energy

rbo(s)

a ine i = 2tt / d& 6 - Trir A H In 5 . (162)
Jo

If the effective radius R (which controls parton shadowing) is energy-independent, the cross section
increases only logarithmically with rising energy. However, the parameter R itself depends on the colli-
sion energy through a convolution with the parton momentum fractions, R 2 (s) ~ Rq + 4 a f eS In 2 5, with
a' ef £ « 0.11 GeV -2 [145]. Thus, the QGSJET cross section exhibits a faster than In s rise,

^inei ~ 4tt ag ff A H In 2 5 . (163)

In SlBYLL-like models, the transverse density distribution is taken as the Fourier transform of the
proton electric form factor, resulting in an energy-independent exponential (rather than Gaussian) fall-off
of the parton density profile for large b,

yb7T

where K 3 denotes the modified Bessel function of the third kind and /j 2 « 0.71GeV 2 [221]. Thus, (159)
becomes

XsH ~e-^ S AH , (165)

and (155) is satisfied when

A H

b (s) = — In s. (166)

Therefore, for SlBYLL-like models, the growth of the inelastic cross section also saturates the In 2 s
Froissart bound [223], but with a multiplicative constant which is larger than the one in QGS JET-like
models

A 2

C^inel ~ 7TC —f In 2 S , (167)

where the coefficient c « 2.5 is found numerically [145].

The main characteristics of the pp cascade spectrum resulting from these choices are readily pre-
dictable: the harder form of the SIBYLL form factor allows a greater retention of energy by the leading
particle, and hence less available for the ensuing shower. Consequently, on average SlBYLL-like models
predict a smaller multiplicity than QGS JET-like models [42].

There are three event generators, sibyll [221], qgsjet [222], and dpmjet [224] which are tai-
lored specifically for simulation of hadronic interactions up to the highest cosmic ray energies. The latest
versions of these packages are sibyll 2.1 [225], qgsjet 11-03 [226], and dpmjet III [227]; respec-
tively. In QGSJET, both the soft and hard processes are formulated in terms of Pomeron exchanges. To

56

describe the minijets, the soft Pomeron mutates into a "semihard Pomeron", an ordinary soft Pomeron
with the middle piece replaced by a QCD parton ladder, as sketched in the previous paragraph. This
is generally referred to as the "quasi-eikonal" model. In contrast, SIBYLL and DPMJET follow a "two
channel" eikonal model, where the soft and the semi-hard regimes are demarcated by a sharp cut in the
transverse momentum: SIBYLL 2.1 uses a cutoff parametrization inspired in the double leading logarith-
mic approximation of the DGLAP equations,

p™ in (^) = p ° T + 0.065 GeV exp[0.9 Vh^} , (168)

whereas dpmjet III uses an ad hoc parametrization for the transverse momentum cutoff

^ in (^) = P° T + 0.12 GeV [log 10 (^/50GeV)] 3 , (169)

where p° T = 2.5 GeV [208].

The transition process from asymptotically free partons to colour-neutral hadrons is described in
all codes by string fragmentation models [228]. Different choices of fragmentation functions can lead
to some differences in the hadron multiplicities. However, the main difference in the predictions of
QGSJET-like and SlBYLL-like models arises from different assumptions in extrapolation of the parton
distribution function to low energy.

The proton-air collisions of interest for air shower development cause additional headaches for
event generators. Both SIBYLL and QGSJET adopt the Glauber formalism [216], which is equivalent to
the eikonal approximation in nucleon-nucleon scattering, except that the nucleon density functions of the
target nucleus are folded with that of the nucleon. The inelastic and production cross sections read:

<ef ir ~ dbb i 1 - ex P K>t AT N (b)]} , (170)

<"od r ~^Jdbb{l- exp [a inel AT N (b)}} , (171)

where T/v(6) is the transverse distribution function of a nucleon inside a nucleus. Here, cr ine i and <7 to t
are given by Eqs. (157) and (158), respectively. The p-air inelastic cross section is the sum of the "quasi-
elastic" cross section, which corresponds to cases where the target nucleus breaks up without production
of any new particles, and the production cross section, in which at least one new particle is generated.
Clearly the development of EASs is mainly sensitive to the production cross section. Overall, the geo-
metrically large size of nitrogen and oxygen nuclei dominates the inclusive proton-target cross section,
and as a result the disagreement from model-dependent extrapolation is not more than about 15%. More
complex nucleus-nucleus interactions are discussed in [229].

EXERCISE 3.1 Consider a typical air nuclei of average (A) = 14.5 and calculate the proton-air
cross section using the approximated expressions (163) and (167) together with the z-integrated Woods-
Saxon profile [230]

T N (b) = | J°° dz {l + exp [(Vb 2 + z 2 - R N )/a\ } 1 ,

where

47T 3

z y N

(172)

(173)

a — 0.5 fm and Rn — 1.1A 1 / 3 fm [231]. Compare the results with those shown in Fig. 25.

57

Table 5: Different k- values used in cosmic ray experiments.

Experiment

k

Fly's Eye

1.6

Akeno

1.5

Yakutsk-99

1.4

EAS-TOP

1.15

Adding a greater challenge to the determination of the proton air cross section at ultrahigh en-
ergies is the lack of direct measuremnts in a controlled laboratory environment. The measured shower
attenuation length, A m , is not only sensitive to the interaction length of the protons in the atmosphere,

Ap-air, With

Am = fcA p _ air = fc^^, (174)

°prod

but also depends on the rate at which the energy of the primary proton is dissipated into EM shower
energy observed in the experiment. Here, A m and A p _ a i r are in gem -2 , the proton mass m p is in g, and
the inelastic production cross section cTp~^ ir is in mb. The value of k depends critically on the inclusive
particle production cross section and its energy dependence in nucleon and meson interactions on the
light nuclear target of the atmosphere. The measured depth X max at which a shower reaches maximum
development in the atmosphere has been the basis of cross section measurements from experiments prior
to HiRes and Auger. However, A max is a combined measure of the depth of the first interaction (which
is determined by the inelastic cross section) and of the subsequent shower development (which has to be
corrected for). The model dependent rate of shower development and its fluctuations are the origin of
the deviation of k from unity in Eq. (174). As can be seen in Table 5, there is a large range of k values
(from 1.6 for a very old model where the inclusive cross section exhibited Feynman scaling, to 1.15 for
modern models with large scaling violations) that make the published values of cr p _ a i r unreliable.

Recently, the HiRes Collaboration developed a quasi-model-free method of measuring o^^ 1 di-
rectly [232]. This is accomplished by folding a randomly generated exponential distribution of first
interaction points into the shower development program, and therefore fitting the entire distribution and
not just the trailing edge. Interestingly, the measured k = 1.21^09 by tne HiRes group is in agreement
with the one obtained by tuning the data to the theory [233, 234].

A compilation of published proton-air cross section measurements is shown in Fig. 25. In the
left panel we show the data without any modification. In the right panel, the published values of o^~^ T
for Fly's Eye [242], Akeno [238], Yakutsk-99 [240], and EAS-TOP [236] collaborations have been
renormalized using the common value of k = 1.264 ± 0.033(stat) ± 0.013(syst) [234]. We have
parametrized the rise of the cross section using a functional form that saturates the Froissart bound,

a P; Q f = A- £ln(£/GeV) + Cln 2 (£/GeV) mb . (175)

The curve with a fast rise, hereafter referred to as case-i, corresponds to A = 280, 5 = 5.7, and C = 0.9.
The slow rise of case-iz has the following parameters: A = 290, B = 6.2, and C = 0.64. The behavior
of the cross section in case-i roughly matches the one implemented in SIBYLL-2.1 [225].

Some guidance towards understanding hadronic processes in the forward direction may come
directly from measurements of hadrons in airshowers [243]. However, the most useful experimental
input in the forseeable future will likely come from the LHC. Processes with low momentum transfer
tend to populate the region at very small angles $ with respect to the beam direction. The distribution
of pseudorapidity, 77 = — lntan($/2), and the energy flow distribution are shown in Fig. 26. While the
particle multiplicity is greatest in the low \rj\ region, it is clearly seen that the energy flow is peaked at

58

Equivalent Vs (GeV) for pp scattering

Equivalent Vs (GeV) for pp scattering

£700

600

500

_0 ARGO-YBJ

1 EAS-TOP
_ A HiRes

A KASCADE prototype
® Akeno
O Yakutsk-90
_ □ Yakutsk-99

O Tien-Shon EAS complex
¥ Fly's Eye

Tevatron LHC

E (GeV)

i0

10

1(T

b700 -

00

_0 ARGO-YBJ

I EAS-TOP (renormalized) Tevatr on
A HiRes

A KASCADE prototype
# Akeno (renormalized)
_0 Yakutsk-90
□ Yakutsk-99 (renormalized)
O Tien— Shan EAS complex
W Fly*s Eye (renormalized)

E (GeV)

Fig. 25: Compilation of proton-air production cross section from cosmic ray measurements (ARGO-YBJ [235],
EAS-TOP [236], HiRes [232], KASCADE prototype [237], Akeno [238], Yakutsk-90 [239], Yakutsk-99 [240],
Tien-Shan EAS complex [241] and Fly's Eye [242]). The data are compared to the parametrizations discussed in
the text; case-i corresponds to the dashed line and case-n to the dot-dashed line.

small production angles (large \rj\). The two general-purpose experiments, ATLAS and CMS, cover up to
\rj\ < 5. TOTEM's coverage in the very forward rapidity range, 3.1 < 77 < 6.5, is committed to measure
the pp total and elastic scattering cross sections, see Appendix C. The fragmentation region that plays
a crucial role in the development of EASs tends to populate the region corresponding to pseudorapidity
range 6 < \rj\ < 10.

The first batch of data collected by the ATLAS experiment has been used to study the rise of the
pp cross section [245]. In Fig. 27 we show a comprehensive investigation into the physics underpin-
ning the rise of cr to t and cr- ine \. 6 The figure shows several models which employ different assumptions
about the evolution of the PDFs at low-x (GRV [249], GRV94 [250], GRV98 [251], MRST [252], and
CTEQ [253]). Model I is an eikonal minijet model incorporating effects induced by the transverse
momentum distribution of soft gluons [254]. In the spirit of the Bloch-Nordsieck study in electrody-
namics [255], the parton distribution in b- space is determined as the Fourier transform of the resummed
distribution of soft gluons (down to zero gluon momenta) emitted from the initial state during the colli-
sion. The very large b limit of the profile function

A(b,s) ~e-( A6 ) 2/ \ (176)

yields an impact parameter distribution falling at its faster like a Gaussian ((3 = 1) and at its slowest like
an exponential ((3 = 1/2). The scale A oc Aqcd includes a mild energy dependence as well as a residual
dependence upon the parameter (3. This behavior in impact parameter space together with the high energy
behavior of the minijet cross section leads to an asymptotic rise of the total cross section [256]

* tot ~f 2 A^ InV^, (177)

consistent with the Froissart bound. In Model II the effective density overlaping parton distributions in
the colliding protons is a SlBYLL-like exponentially decreasing hyperbolic Bessel function, determined

6 The CMS measurements of the pp cross section at y/s = 7 TeV (not shown in the figure) are summarized in [248].

59

I§" 8
O

? 7

6
5
4
3
2
1

W 1.4

1 .2

0.8
0.6
0.4
0.2

I i i i i I

<

LHC, inelastic collisions

TOTEM, CMS

Tracking and Calorimetry

ATLAS, CMS

< 7-7—, * ; >

>

Calorimetry

ATLAS, CMS
Tracking

iii

j— i i i i i — i — \ i i i i i i i i i i i i i i i i i i

-2.5 -5 -2.5 2.5 5 2.5 10

Fig. 26: Pseudorapidity distributions of charged particles (upper panel) and of the energy flow (lower panel) for
pp collisions at LHC. From Ref. [244].

from the proton electric form factor [257]. The analyticity-constrained analytic amplitude model of
Block and Halzen [258, 259] exploits high quality low energy data to determine an asymptotic growth
of the cross section, <7 to t ~ In 2 s, which saturates the Froissart bound. The Aspen model is a revised
version of the eikonal model of Block, Gregores, Halzen, and Pancheri [260], which now incorporates
analyticity constraints [261].

In a complementary study, a comparison of LHC data (on inclusive particle production at ener-
gies = 0.9, 2.36, and 7 TeV) with predictions of a variety of hadronic interaction models has been
elaborated in [262]. The results show that while reasonable qualitative agreement has been achieved for
some event generators, none of them reproduces the evolution of all the observables with particularly
compeling precision.

In summary, high energy hadronic interaction models are still being refined and therefore the
disparity between them can vary even from version to version. At the end of the day, however, the relevant
parameters boil down to two: the mean free path, AcR-air = (^p r od n atm) _1 , and the inelasticity,
^cR-air = 1 — ^lead/^proj. where n atm is the number density of atmospheric target nucleons, £i e ad is the
energy of the most energetic hadron with a long lifetime, and E pro j is the energy of the projectile particle.
The first parameter characterizes the frequency of interactions, whereas the second one quantifies the
energy lost per collision. Overall, SIBYLL has a shorter mean free path and a smaller inelasticity than

60

E

100

120

80

a tot Model I
a inel Model II

_ _ _ GRV, pTmin=1 .15 GeV, a T0T (7 TeV)=93.4 mb, a INEL (7 TeV)=56.3 mb

* Aspen Model & Analytic Amplitude Model

* ATLAS - PYTHIA extrapolation - preliminary

* ATLAS - preliminary data

i

60

20

40

• UA5

A UA1

O UA4

* CDF
□ E710
▼ E811

4

Vs ( GeV )

10

Fig. 27: Compilation of total and inelastic pp and pp cross section measurements. The data are compared to model
predictions. The dashed lines indicate the predictions of Model II [246]. From low to high energy the same eikonal
function has been used to compute both the total and the inelastic cross sections. The predictions are calculated
using GRV pfd's, with p™ in = 1.15 GeV. The green band is obtained with Model I varying by no more
than a few percent and using different partonic densities [246]. The blue band corresponds to Model II using
pf in = 1.1 GeV (with MRST) in the upper edge and p™ in = 1.15 GeV (with GRV94) in the lower edge [246].
Predictions for a tot from the analytic amplitude model and for cr inel from the Aspen model are indicated by blue
stars [247]. This figure is courtesy of Giulia Pancheri.

QGSJET. Since a shorter mean free path tends to compensate for a smaller inelasticity, the two codes
generate similar predictions for an air shower which has lived through several generations. Both models
predict the same multiplicity below about 10 7 GeV, but the predictions diverge above that energy. Such
a divergence readily increases with rising energy. While QGSJET predicts a power law-like increase of
the number of secondaries up to the highest energy, SIBYLL multiplicity exhibits a logarithmic growth.
As it is extremely difficult to observe the first interactions experimentally, it is not straightforward to
determine which model is closer to reality.

3.2 Electromagnetic processes

The evolution of an extensive air shower is dominated by electromagnetic processes. The interaction of
a baryonic cosmic ray with an air nucleus high in the atmosphere leads to a cascade of secondary mesons
and nucleons. The first few generations of charged pions interact again, producing a hadronic core, which
continues to feed the electromagnetic and muonic components of the showers. Up to about 50 km above
sea level, the density of atmospheric target nucleons is n atm ~ 10 20 cm -3 , and so even for relatively
low energies, say E n ± « 1 TeV, the probability of decay before interaction falls below 10%. Ultimately,
the electromagnetic cascade dissipates around 90% of the primary particle's energy, and hence the total
number of electromagnetic particles is very nearly proportional to the shower energy [263].

Roughly speaking, at 10 11 GeV, baryons and charged pions have interaction lengths of the order of
40 g/cm 2 , increasing to about 60 g/cm 2 at 10 7 GeV. Additionally, below 10 10 GeV, photons, electrons,
and positrons have mean interaction lengths of 37 g/cm 2 . Altogether, the atmosphere acts as a natural

61

colorimeter with variable density, providing a vertical thickness of 26 radiation lengths and about 15
interaction lengths. Amusingly, this is not too different from the number of radiation and interaction
lengths at the LHC detectors. For example, the CMS electromagnetic calorimeter is > 25 radiation
lengths deep, and the hadron calorimeter constitutes 1 1 interaction lengths.

By the time a vertically incident 10 11 GeV proton shower reaches the ground, there are about 10 11
secondaries with energy above 90 keV in the the annular region extending 8 m to 8 km from the shower
core. Of these, 99% are photons, electrons, and positrons, with a typical ratio of 7 to e + e~ of 9 to 1.
Their mean energy is around 10 MeV and they transport 85% of the total energy at ground level. Of
course, photon-induced showers are even more dominated by the electromagnetic channel, as the only
significant muon generation mechanism in this case is the decay of charged pions and kaons produced in
7-air interactions [264].

It is worth mentioning that these figures dramatically change for the case of very inclined showers.
For a primary zenith angle, 9 > 70° , the electromagnetic component becomes attenuated exponentially
with atmospheric depth, being almost completely absorbed at ground level. Note that the vertical atmo-
sphere is « 1, 000 g/cm 2 , and is about 36 times deeper for completely horizontal showers. As a result,
most of the energy at ground level from an inclined shower is carried by muons.

In contrast to hadronic collisions, the electromagnetic interactions of shower particles can be cal-
culated very accurately from quantum electrodynamics. Electromagnetic interactions are thus not a
major source of systematic errors in shower simulations. The first comprehensive treatment of electro-
magnetic showers was elaborated by Rossi and Greisen [265]. This treatment was recently cast in a more
pedagogical form by Gaisser [63], which we summarize in the subsequent paragraphs.

The generation of the electromagnetic component is driven by electron bremsstrahlung and pair
production [266]. Eventually the average energy per particle drops below a critical energy, eo, at which
point ionization takes over from bremsstrahlung and pair production as the dominant energy loss mech-
anism. The e ± energy loss rate due to bremsstrahlung radiation is nearly proportional to their energy,
whereas the ionization loss rate varies only logarithmically with the e ± energy. Though several differ-
ent definitions of the critical energy appear in the literature [267], throughout these lectures we take the
critical energy to be that at which the ionization loss per radiation lenght is equal to the electron energy,
yielding e = 710 MeV/(Z eff + 0.92) ~ 86 MeV [268]. 7 The changeover from radiation losses to ion-
ization losses depopulates the shower. One can thus categorize the shower development in three phases:
the growth phase, in which all the particles have energy > eo; the shower maximum, X max ; and the
shower tail, where the particles only lose energy, get absorbed or decay.

The relevant quantities participating in the development of the electromagnetic cascade are the
probability for an electron of energy E to radiate a photon of energy k = y hrem E and the probability
for a photon to produce a pair e + e~ in which one of the particles (hereafter e~) has energy E = y pair k.
These probabilities are determined by the properties of the air and the cross sections of the two processes.

In the energy range of interest, the impact parameter of the electron or photon is larger than an
atomic radius, so the nuclear field is screened by its electron cloud. In the case of complete screen-
ing, where the momentum transfer is small, the cross section for bremsstrahlung can be approximated

where A e $ is the effective mass number of the air, X EM is a constant, and is Avogadro's number. In
the infrared limit (i.e. y hrem <C 1) this approximation is inaccurate at the level of about 2.5%, which is
small compared to typical experimental errors associated with cosmic air shower detectors. Of course,

7 For altitudes up to 90 km above sea level, the air is a mixture of 78.09% of N2, 20.95% of O2, and 0.96% of other
gases [269]. Such a mixture is generally modeled as an homogeneous substance with atomic charge and mass numbers Z e ff =
7.3 and A e s = 14.6, respectively.

by [270]

(178)

62

the approximation fails as y hvem — » 1, when nuclear screening becomes incomplete, and as y hrem — >> 0, at
which point the LPM and dielectric suppression effects become important, as we discuss below.

Using similar approximations, the cross section for pair production can be written as [270]

<fo 7 -> e + e - _ Aeff ( 4 4 2 \

d# ~ x em tv a v 3?/pair 3 ^ pair /'

The similarities between this expression and Eq. (178) are to be expected, as the Feynman diagrams for
pair production and bremsstrahlung are variants of one another.

The probability for an electron to radiate a photon with energy in the range (fc, k+dk) in traversing
dt = dX/X BM of atmosphere is

da e ^ X EM N A A+ _( ni 41- 2/ brem

whereas the corresponding probability density for a photon producing a pair, with electron energy in the
range (E, E + dE), is

da 7 ^ e+e - X EM 7V A / 4 4 2

« 1 - -y . + . ) dy . dt . (181)

dE A eff V 3 3 pai V

The total probability for pair production per unit of X EM follows from integration of Eq. (181),

7 d# A eff y V 3 pair 3^ air y ^ pair 9 v ;

As can be seen from Eq. (180), the total probability for bremsstrahlung radiation is logarithmically
divergent. However, this infrared divergence is eliminated by the interference of bremsstrahlung ampli-
tudes from multiple scattering centers. This collective effect of the electric potential of several atoms
is known as the LPM effect [32, 33]. Of course, the LPM suppression of the cross section results in an
effective increase of the mean free path of electrons and photons. This effectively retards the develop-
ment of the electromagnetic component of the shower. It is natural to introduce an energy scale, ^lpm ?
at which the inelasticity is low enough that the LPM effect becomes significant [271]. Below £?lpm, the
energy loss rate due to bremsstrahlung is roughly

dE
dX

f 2/ brem E L m + \ 1 J hre A dy hrem =-^~. (183)

A EM ^0 V 6 Vbrem J A EM

With this in mind, we now identify the constant X EM « 36.7 g cm - ^ with the radiation length in air
defined as the mean distance over which a high-energy electron loses 1 je of its energy, or equivalently
7/9 of the mean free path for pair production by a high-energy photon [267].

The most evident signatures of the LPM effect on shower development are a shift in the position
of the shower maximum X max and larger fluctuations in the shower development. When considering
the LPM effect in the development of air showers produced by UHECRs, one has to keep in mind that
the suppression in the cross sections is strongly dependent on the atmospheric depth. 8 Since the upper
atmosphere is very thin, the LPM effect becomes noticeable only for photons and electrons with energies
above £^lpm ~ 10 10 GeV. For baryonic primaries the LPM effect does not become important until the
primary energy exceeds 10 12 GeV. This is because the electromagnetic shower does not commence until
after a significant fraction of the primary energy has been dissipated through hadronic interactions. To
give a visual impression of how the LPM effect slows down the initial growth of high energy photon-
induced showers, we show the average longitudinal shower development of 10 10 GeV proton and 7-ray
showers (generated using AIRES 2.6.0 [273]) with and without the LPM effect in Fig. 28.

8 The same occurs for dielectric suppression, although the influence is not as important as for the LPM effect [272].

63

x 10
| 3000

2500

2000

1500

1000

500

200 400 600 800 1000 1200 1400 1600 1800

slant depth (gxm =2 )

Fig. 28: Average longitudinal shower developments of 10 11 GeV proton (dashed-dotted line) and 7-rays with and
without the LPM effect (solid and dotted lines, respectively). The primary zenith angle was set to 6 = 60°. The
shadow area represents the intrinsic fluctuations of the showers. Larger fluctuations can be observed for 7-ray
showers with the LPM effect, as expected. From Ref . [200] .

At energies at which the LPM effect is important {viz. E > ElpmX 7-ray showers will have al-
ready commenced in the geomagnetic field at almost all latitudes: primary photons with E > 10 10 GeV
convert into e + e~ pairs, which in turn emit synchrotron photons. This reduces the energies of the pri-
maries that reach the atmosphere, and thereby compensates for the tendency of the LPM effect to retard
the shower development. The relevant parameter to determine both conversion probability and syn-
chrotron emission is E x B±, where E is the 7-ray energy and B± the transverse magnetic field. This
leads to a large directional and geographical dependence of shower observables. Thus, each experiment
has its own preferred direction for identifying primary gamma rays. For instance, Fig. 29 shows a map
of the photon conversion probability in the geomagnetic field for all incident directions evaluated at the
location of the HiRes experiment (\B\ = 0.53 G, 1 = 25°, and S = 14°) [274]. The smallest probabilities
for conversion are found, not surprisingly, around the direction parallel to the local geomagnetic field.
Note that this conversion-free region shrinks rapidly with increasing primary energy. A similar evalua-
tion for the Southern Site of the Pierre Auger Observatory (\B\ = 0.25 G, 1 = —35°, and S = 86°) can
be found in [275].

The muonic component of an EAS differs from the electromagnetic component for two main
reasons. First, muons are generated through the decay of cooled (E^± < 1 TeV) charged pions, and thus
the muon content is sensitive to the initial baryonic content of the primary particle. Furthermore, since
there is no "muonic cascade", the number of muons reaching the ground is much smaller than the number
of electrons. Specifically, there are about 5 x 10 8 muons above 10 MeV at ground level for a vertical
10 11 GeV proton induced shower. Second, the muon has a much smaller cross section for radiation and

64

E=10 195 [eV] E=10 200 [eV] E=10 205 [eV] E=10 210 [eV]

50 100% 50 100% o 50 100% 50 100%

Conversion Probability Conversion Probability Conversion Probability Conversion Probability

Fig. 29: Maps of gamma ray conversion probability in the geomagnetic field for several primary energies. Az-
imuths are as labeled, "N" denotes true north. The inner circles correspond to zenith angles 30°, 60° and horizon.
Dashed curves indicate the opening angles of 30°, 60° and 90° to the local magnetic field. From Ref. [274].

pair production than the electron, and so the muonic component of an EAS develops differently than
does the electromagnetic component. The smaller multiple scattering suffered by muons leads to earlier
arrival times at the ground for muons than for the electromagnetic component.

The ratio of electrons to muons depends strongly on the distance from the core; for example,
the e + e~ to ratio for a 10 11 GeV vertical proton shower varies from 17 to 1 at 200 m from the

core to 1 to 1 at 2000 m. The ratio between the electromagnetic and muonic shower components behaves
somewhat differently in the case of inclined showers. For zenith angles greater than 60°, the e + e~ I /x + \±~
ratio remains roughly constant at a given distance from the core. As the zenith angle grows beyond
60°, this ratio decreases, until at 9 = 75°, it is 400 times smaller than for a vertical shower. Another
difference between inclined and vertical showers is that the average muon energy at ground level changes
dramatically. For horizontal showers, the lower energy muons are filtered out by a combination of energy
loss mechanisms and the finite muon lifetime: for vertical showers, the average muon energy is 1 GeV,
while for horizontal showers it is about 2 orders of magnitude greater. The muon densities obtained in
shower simulations using SIBYLL 2.1 fall more rapidly with lateral distance to the shower core than those
obtained using QGSJET 01. This can be understood as a manifestation of the enhanced leading particle
effect in SIBYLL, which can be traced to the relative hardness of the electromagnetic form factor profile
function. The curvature of the distribution ((Pp^/dr 2 ) is measurably different in the two cases, and,
with sufficient statistics, could possibly serve as a discriminator between hadronic interaction models,
provided the primary species can be determined from some independent observable(s) [276].

High energy muons lose energy through pair production, muon-nucleus interaction, bremsstrahlung,
and knock-on electron (5-ray) production [277]. The first three processes are discrete in the sense that
they are characterized by high inelasticity and a large mean free path. On the other hand, because of its
short mean free path and its small inelasticity, knock-on electron production can be considered a con-
tinuous process. The muon bremsstrahlung cross section is suppressed by a factor of {m e /m^) 2 with
respect to electron bremsstrahlung; see Eq. (178). Since the radiation length for air is about 36.7 g/cm 2 ,
and the vertical atmospheric depth is 1, 000 g/cm 2 , muon bremsstrahlung is of negligible importance
for vertical air shower development. Energy loss due to muon-nucleus interactions is somewhat smaller
than muon bremsstrahlung. As can be seen in Fig. 30, energy loss by pair production is slightly more
important than bremsstrahlung at about 1 GeV, and becomes increasingly dominant with energy. Finally,
knock-on electrons have a very small mean free path (see Fig. 30), but also a very small inelasticity, so
that this contribution to the energy loss is comparable to that from the hard processes.

3.3 Paper-and-pencil air shower modeling

Most of the general features of an electromagnetic cascade can be understood in terms of the toy model
due to Heitler [278]. In this model, the shower is imagined to develop exclusively via bremsstrahlung

65

and pair production, each of which results in the conversion of one particle into two (see Fig. 31). As
was previously discussed, these physical processes are characterized by an interaction length X EM w
37.6 g/cm 2 . One can thus imagine the shower as a particle tree with branches that bifurcate every X EM ,
until they fall below a critical energy, eo ~ 86 MeV, at which point energy loss processes dominate.
Up to eo, the number of particles grows geometrically, so that after n = X/X EM branchings, the total
number of particles in the shower is N w 2 n . At the depth of shower maximum X max , all particles are
at the critical energy, eo, and the energy of the primary particle, Eq, is split among all the Af max = Eq/co
particles. Putting this together, we get:

ln(£ /e )
In 2

^max ~ ^em — ttt; — ■ (184)

Changes in the mean mass composition of the CR flux as a function of energy will manifest as changes
in the mean values of X max . This change of X max with energy is commonly known as the elongation
rate [279]:

For purely electromagnetic showers, X max (£?) « X EM \n(E/eo), and hence D e « Xem- For con-
venience, the elongation rate is often written in terms of energy decades, Dio = <9(X max ) / 'dlogE,
where Dio = 2.3D e . The elongation rate obtain from the Heitler model, Dj « 84 g/cm 2 , is in very
good agreement with the results from Monte Carlo simulations. However, the prediction for the particle
number at maximum is overestimated by a factor of about 2 to 3. Moreover, Heitler' s model predicts a
ratio of electron to photons of 2, whereas simulations and direct cascade measurements in the air show

66

n=l

n=2

n=3

(b)

P

n=

n=4

n=2

n=3

Fig. 31: Schematic views of (a) an electromagnetic cascade and (b) a hadronic shower. In the hadron shower,
dashed lines indicate neutral pions which do not re-interact, but quickly decay, yielding electromagnetic subshow-
ers (not shown). Not all pions lines are shown after the n = 2 level. Neither diagram is to scale. From Ref. [280].

a ratio of the order of 1/6. This difference is due to the fact that multiple photons are emitted during
bremsstrahlung and that electrons lose energy much faster than photons do.

As we have seen, baryon-induced showers are also dominated by electromagnetic processes, thus
Heitler's toy model is still enlightening for such cases. For proton primaries, the multiplicity rises with
energy, and the resulting elongation rate becomes smaller. This can be understood by noting that, on
average, the first interaction is determined by the proton mean free path in the atmosphere, A p _ a i r = Xq.
In this first interaction the incoming proton splits into (n(E)) secondary particles, each carrying an
average energy E/{n(E)). Assuming that X max (E) depends dominantly on the first generation of 7
subshowers, the depth of maximum is obtained as in Eq. (184),

X max (£) « X + X EM \u[E/{n(E))]

(186)

For a proper evaluation of X max , it would be necessary to sum each generation of subshowers carefully
from their respective points of origin, accounting for their attenuation near and after the maxima. If we
now further assume a multiplicity dependence (n(E)) « noE A , then the elongation rate becomes,

5X„

SlnE

1

Sln(n{E))
SlnE

+

sx

SlnE

which corresponds to the form given by Linsley and Watson [31],

D e = I c

1 -

Sln{n{E)) X 6\n(X )

SlnE

+

*™ SlnE

= X EM (! _ B ) >

where

B = A —

X SlnX

X^„ SlnE

(187)

(188)

(189)

A precise calculation of a proton shower evolution has been carried out by Matthews [280], using
the simplifying assumption that hadronic interactions produce exclusively pions, 2N W charged and
neutral (see Fig. 31). 7r°'s decay immediately and feed the electromagnetic component of the shower,
whereas ir^ 's soldier on. The hadronic shower continues to grow, feeding the electromagnetic component
at each interaction, until charged pions reach a characterictic energy at which decay is more likely than
a new interaction. The interaction length and the pion multiplicity (3N W ) are energy independent in the
Heitler-Matthews model. The energy is equally shared by the secondary pions. For pion energy between
1 GeV and 10 TeV, a charged multiplicity of 10 (N w = 5) is an appropriate number.

67

The first interaction diverts 1/3 of the available energy (Eq/3) into the EM component while the
remaining 2/3 continue as hadrons. The number of hadrons increases through subsequent generation of
particles and in each generation about 30% of the energy is transferred to the EM cascade. Therefore
the longer it takes for pions to reach the characteristic energy £J ~ 20 GeV (below which they will
decay into muons), the larger will be the EM component. Consequently, in long developing showers the
energy of the muons from decaying pions will be smaller. In addition, because of the density profile of
the atmosphere, ^ is larger high above ground than at sea level and deep showers will produce fewer
muons.

This positive correlation introduces a link between the primary cosmic ray interaction cross section
on air and the muon content at ground level. According to those principles, primaries with higher cross
sections will have a larger muon to electron ratio at ground level.

To obtain the number of muons in the shower, one simply assumes that all pions decay into muons
when they reach the critical energy: = (2N 7T ) Uc , where n c = \u(Eq/^) / ln(3A^-) is the number of
steps needed for the pions to reach £J . Introducing (3 = ln(2A^-) / ln(3A^-) we have

N» = (E /g)P . (190)

For N n = 5, (3 = 0.85. Unlike the electron number, the muon multiplicity does not grow linearly with
the primary energy, but at a slower rate. The precise value of (3 depends on the average pion multiplicity
used. It also depends on the inelasticity of the hadronic interactions. Assuming that only half of the
available energy goes into the pions at each step (rather than all of it, as done above) would lead to
(3 = 0.93. Detailed simulations give values of f3 in the range 0.9 to 0.95 [145].

The first interaction yields N 1 = 2A^.o = N n ± . Each photon initiates an EM shower of energy
Eo/(3N 7r ±) = Eq/(6N 7T ). Using pp data [15], we parametrized the charged particle production in the
first interaction as N n ± = 41.2(£'o/l PeV) 1 / 5 . Now, from the approximation in (186), based on the sole
evolution of the EM cascade initiated by the first interaction, we obtain

XL* = Xo + X^HEo/iGN^eo)}

= (470 + 581ogio[£;o/lPeV])g/cm 2 . (191)

This falls short of the full simulation value by about 100 g/cm 2 [280].

EXERCISE 3.2 The depth of shower maximum obtained in Eq. (191) is only approximate since
it considers just the first interaction as hadronic in nature. Extend the approximation to include hadronic
interactions in the second generation of particles. Try the generalization to include all generations of
hadronic collisions until the charged pions cool down below the critical energy.

A good approximation of the elongation rate can be obtained when introducing the cross section
and multiplicity energy dependence. Using a p-air cross section of 550 mb at 10 9 GeV and a rate of
change of about 50 mb per decade of energy leads to [281]

X ~ 90 - 9 log (Eo/EeV) g/cm 2 . (192)

Now, assuming (as in [280]) that the first interaction initiates 2N n EM cascades of energy Eq/QN^, with
N n oc (^o/PeV) 1 / 5 for the evolution of the first interaction multiplicity with energy, we can calculate
the elongation rate

yp ^ ^max = d(X In2 + X EM ln[£ /(6A^e )]
10 dlogE dlog^o

jjP — -iiidA -y— u — ~ ■ -EM L u / V— (193)

or

5

Dl Q = ^Dj - 9 In 2 - 62 g/cm 2 . (194)

68

This result is quite robust as it only depends on the cross section and multiplicity evolution with energy.
It is in good agreement with Monte Carlo simulation [145].

To extend this discussion to heavy nuclei, we can apply the superposition principle as a reasonable
first approximation. In this approximation, we pretend that the nucleus comprises unbound nucleons,
such that the point of first interaction of one nucleon is independent of all the others. Specifically, a
shower produced by a nucleus with energy E A and mass A is modeled by a collection of A proton
showers, each with A -1 of the nucleus energy. Modifying Eq. (184) accordingly one easily obtains
I max oc \u(Eq/A). Assuming that B is not changing with energy, one obtains for mixed primary
composition [31]

d(]nA)'

D e = X (1 - B)

1 -

d\nE

(195)

Thus, the elongation rate provides a measurement of the change of the mean logarithmic mass with
energy. One caveat of the procedure discussed above is that Eq. (188) accounts for the energy dependence
of the cross section and violation of Feynman scaling only for the first interaction. Note that subsequent
interactions are assumed to be characterized by Feynman scaling and constant interaction cross sections;
see Eq. (189). Above 10 7 GeV, these secondary interactions play a more important role, and thus the
predictions of Eq. (195), depending on the hadronic interaction model assumed, may vary by up to
20% [145]. In Fig. 6 we show the variation of (X max ) with primary energy as measured by HiRes and
Auger together with predictions from a variety of hadronic interaction models.

The muon content of an EAS at ground level 7V M , as well as the ratio N^/N e , are sensitive to
primary composition (here, N e is the electron content at ground level). To estimate the ratio of the
muon content of nucleus-induced to proton-induced showers, we can resort again to the principle of
superposition. Using f3 = 0.93 we find that the total number of muons produced by the superposition of
A individual proton showers is, oc A(E A /A) - 93 . Consequently, in a vertical shower, one expects
a cosmic ray nucleus to produce about A om more muons than a proton. This implies that a shower
initiated by an iron nucleus produces about 30% more muons than a proton shower. Note, however,
that a change in the hadronic interaction model could produce a much larger effect than a change in the
primary species. For example, replacing QGSJET 01 with SIBYLL 1.6 as the hadronic interaction model
leads to a prediction of 60% more muons for an iron shower than for a proton shower [282].

The situation for gamma-induced showers is a bit different. In this case the muon component of
the shower does not simply follow Eq. (190) because of the LPM and geomagnetic field effects [283].
Competition between the two processes leads to a complex behavior in N^/Nfl, as shown in Fig. 32.

While these toys models are very useful for imparting a first intuition regarding global shower
properties, the details of shower evolution are far too complex to be fully described by a simple analytical
model. Full Monte Carlo simulation of interaction and transport of each individual particle is required
for precise modeling of the shower development. At present, two Monte Carlo packages are available
to simulate EASs: CORSIKA (COsmic Ray Simulation for KAscade) [284] and AIRES (AIR shower
Extended Simulation) [273]. Both programs provide fully 4-dimensional simulations of the air showers
initiated by protons, photons, and nuclei. To simulate hadronic physics, the programs make use of the
event generators described in Sec. 3.1. Propagation of particles takes into account the Earth's curvature
and geomagnetic field. For further details on these codes, the reader is referred to [285].

4 Searches for new physics beyond the electro weak scale at a/s ^ 250 TeV
4.1 General idea

If new physics interations occur at LHC energies, then CR collisions with cm. energies ranging up to
250 TeV would obviously involve new physics as well. The question is, can new physics be detected by
CR experiments? At ultrahigh energies, the cosmic ray luminosity ~7x 10 -10 (E/PeV)~ 2 cm -2 s _1
(taking a single nucleon in the atmosphere as a target and integrating over 2tt sr) is about 50 orders of

69

0.04

i i 1 1 1 1 ill

I I I M ill I I I I M ill I I I I M ill I I I llllll

I I llllll I I I I

10

17

10

10

19

10

20

10

21

10

22

E(eV)

Fig. 32: Ratio of the muon content for EASs produced by primary gammas and protons. The geomagnetic field is
set to the PAO Southern site. From Ref. [283].

magnitude smaller than the LHC luminosity. This renders the hunt for physics beyond the electroweak
scale futile in hadronic cosmic ray interactions occuring in the atmosphere. However, there is still a
possibility of uncovering new physics at sub-fermi distances in cosmic neutrino interactions.

Neutrinos are unique probes of new physics, as their interactions are uncluttered by the strong and
electromagnetic forces and, upon arrival at the Earth, they may experience collisions with cm. energies
up to < 250 TeV. However, rates for new physics processes are difficult to test since the flux of cosmic
neutrinos is virtually unknown. Interestingly, it is possible in principle to disentangle the unknown flux
and new physics processes by using multiple observables [286,287].

For example, possible deviations of the neutrino-nucleon cross section due to new non-perturbative
interactions 9 can be uncovered at Auger by combining information from Earth- skimming and quasi-
horizontal showers. In particular, if an anomalously large rate is found for deeply developing quasi-
horizontal showers, it may be ascribed either to an enhancement of the incoming neutrino flux, or
an enhancement in the neutrino-nucleon cross section (assuming non-neutrino final states dominate).
However, these possibilities can be distinguished by comparing the rates of Earth- skimming and quasi-
horizontal events. For instance, an enhanced flux will increase both quasi-horizontal and Earth- skimming

9 Throughout these lectures we use this term to describe neutrino interactions in which the final state energy is dominated
by the hadronic component. We are not considering here new "perturbative" physics e.g. (softly broken) supersymmetry at the
TeV scale which would have quite different signatures in cosmic ray showers.

70

event rates, whereas an enhanced interaction cross-section will also increase the former but suppress the
latter, because the hadronic decay products cannot escape the Earth's crust. Essentially this approach
constitutes a straightforward counting experiment, as the detailed shower properties are not employed to
search for the hypothesized new physics.

The question we would like to answer is then how many Earth- skimming and quasi-horizontal
events would we need to observe at Auger to make a convincing case for the existence of non-perturbative
physics in which the final state is dominated by hadrons. The analysis techniques described herein consti-
tute an entirely general approach to searching for non-perturbative interactions without any dependence
on what hypothetical mechanism might actually cause the "hadrophilia." In Appendix D we illustrate
one possible new physics process which may be accessible using these techniques at Auger, and which
is beyond the reach of the LHC.

4.2 v acceptance and systematic uncertainties

Detailed Monte Carlo simulations are used to compute the acceptance for Earth- skimming (ES) and
quasi-horizontal (QH) events. Neutrinos are propagated through the atmosphere, the Earth's crust, and
the Andes mountains using an extended version [288] of the code ANIS [289]. In the simulations, the
vN cross-sections from reference [290] are employed. Particles resulting from vN interactions are fed
to PYTHIA [291] and r decays are simulated using TAUOLA [292].

The flux, energy and decay vertices of outgoing leptons are calculated inside an "active detector"
volume of 3, 000 x 10 km 3 , including the real shape of the surface array. A relief map of the Andes
mountains was constructed using digital elevation data from the Consortium for Spatial Information
(CGIAR-CSI) [293]. The map of the area around the Auger site is depicted in Fig. 33.

Fig. 33: Topography in the vicinity of the Auger site. The surface array is centered at X = Y = 0. From Ref. [294].

To study the response of the detector, the outputs of PYTHIA and/or TAUOLA are used as input
for the AIRES [295] air shower simulation package. The response of the surface detector array is sim-
ulated in detail using the Auger Offline simulation package [296]. Atmospheric background muons are
also simulated in order to study the impact on neutrino identification, as such accidental muons can be
wrongly classified as shower particles. The background from hadronic showers above 10 8 GeV is esti-
mated to be 0(1) in 20 years, so for the energy bin considered in this analysis, 9.5 < log 10 (^/GeV),
the background is negligible.

71

To establish benchmark neutrino rates, we use the Waxman-Bahcall bound [144] for the flux,
<3>wb = 2.3 x 10 -8 E~ 2 GeV -1 s _1 cm -2 sr _1 , and employ the acceptance computed by the simulations
described above. In order to estimate the systematic uncertainty associated with our lack of knowledge
of the dependence of the flux on energy, we consider several scenarios which plausibly bracket the range
of possibilities:

(C/E )E-\
C Eq e~ 3 ,

CE~ 2 exp[- log 10 (E u /E ) 2 /(2a 2 )},

where C = 2.3 x 1CT 8 GeV -1 s -1 cm- 2 sr" 1 , E = 10 10 GeV, a = 0.5 GeV. The expected rates for the
entire range over which Auger is sensitive are given in Table 6 and the rates for the high energy bin are
given in Table 7.

Table 6: Expected events per year at Auger in the energy range 8 < log 10 (£^/GeV), for various incident
zenith angle (0) ranges, assuming the Waxman-Bahcall flux.

flux

up-going

down-going

ratio

(2)
(2)

9 N Ut

9 N Ue N Ut N V/l N VM

90-95 0.68
90-95 0.68

60-90 0.134 0.109 0.019 0.262
75-90 0.075 0.071 0.011 0.157

2.58
4.27

Table 7: Expected events per year (Ni) at Auger in the energy range 9.5 < log 10 (^/GeV) < 10.5, for various
incident zenith angle (0) ranges and the 4 flux models considered.

flux

up-going

down-going

ratio

9

9

(1)

90-95

0.14

60-90

0.059

0.049

0.011

0.12

1.14

(2)

90-95

0.15

60-90

0.059

0.049

0.096

0.11

1.33

(3)

90-95

0.23

60-90

0.079

0.062

0.0123

0.15

1.53

(4)

90-95

0.12

60-90

0.046

0.037

0.0080

0.091

1.33

(1)

90-95

0.14

75-90

0.027

0.031

0.0056

0.064

2.14

(2)

90-95

0.15

75-90

0.026

0.029

0.0048

0.060

2.47

(3)

90-95

0.23

75-90

0.036

0.041

0.0062

0.083

2.75

(4)

90-95

0.12

75-90

0.021

0.024

0.0040

0.049

2.45

Table 8 contains a summary of systematic uncertainties on the ratio of the number of ES to QH
events. The uncertainty in spectrum shape is taken from Table 7. The uncertainty on the parton structure
of the nucleon is estimated by considering different PDFs (GRV92NLO [297] and CTEQ66c [298]).
Finally, the uncertainty on the energy loss, /3 r , of r leptons as they propagate through the Earth's crust is
derived from [299-302].

4.3 Auger discovery reach

Consider first a flux of Earth- skimming r neutrinos with energy in the range 10 9 5 GeV < E v <
10 10 - 5 GeV. The neutrinos can convert to r leptons in the Earth via the charged current interaction

1. ^ B (E U ) =

2. $% B (E U ) =

3. ^ B (E U ) =

4. ^ B (E U ) =

72

Table 8: Contributions to the systematic uncertainty on the Earth- skimming to quasi-horizontal event ratio. We
have considered the energy range 9.5 < log 10 (E u / GeV) < 10.5 and the zenith angle range 75° < 6 < 90°.

ratio

flux

PDF

sum

+11%

0%

+24%

+ 26%

2.47

2.47

-13%

-21%

-25%

-35%

v T ± N — >► r^X. In the (perturbative) SM, the interaction path length for the neutrino is

L£ c = [NAPsa^c}- 1 , (196)

where a^ c is the charged current cross section for a neutrino energy E v = Eq. The density of the
material through which the neutrinos pass, p s , is about 2.65 g/cm 3 for the Earth's crust. Here we
have neglected neutral current interactions, which at these energies only reduce the neutrino energy
by approximately 20%, which is within the systematic uncertainty. For Eo ~ 10 10 GeV, L£ c ~
(9(100) km. Let us assume some hypothetical non-perturbative physics process enhances the vN cross
section. Then the interaction path length becomes

L t l ' ot = [iV A p s (^ c + ^p)]- 1 , (197)

where a^ p is the non-perturbative contribution to the cross section for E v = Eq.

Once a r is produced by a CC interaction, it can be absorbed in the Earth or escape and possibly
decay, generating a detectable air shower. At these high energies, the r propagation length in the Earth
is dominated by energy loss rather than the finite r lifetime. The energy loss can be expressed as

^ = -(a T + p T E T )p 8 , (198)

where a T characterizes energy loss due to ionization and f3 T characterizes losses through bremsstrahlung,
pair production and hadronic interactions. At these energies, energy losses due to ionization turn out to
be negligible, while f3 T ~ 0.8 x 10 -6 cm 2 /g [303]. From Eq. (198), we observe that the maximum path
length for a detectable r can be written

L T = ^ln(£ max /£ min ) , (199)

PtPs

where £ max « Eo is the energy at which the r is created, and i? m in is the minimal energy at which a r
can produce a shower big enough to be detected. For ^max/^min = 10, L T = 11 km.

The probability for a neutrino with incident nadir angle 9 to emerge as a detectable r is

pi j 7

P W = / T^ e ~ Z ' LXot 9 [* " ( £ ~ LT )1 ' ( 20 °)
Jo L cc

where £ = 2i? e cos 9 is the chord length of the intersection of the neutrino's trajectory with the Earth,
with i?0 « 6371 km the Earth's radius. Note that we have neglected the possibility that non-perturbative
processes could lead to a detectable signal, since the hadrons which dominate the final state will be
absorbed in the Earth. The step function in Eq. (200) reflects the fact that a r will only escape the Earth
if z + L T > £, as illustrated in Fig 34.

73

Fig. 34: The chord length of the intersection of a neutrino with the Earth is i = 2R§ cos 0. In the figure, the
neutrino produces a lepton / after traveling some distance z inside the Earth's crust. If z + L r > £, the lepton will
escape the Earth and can generate an air shower. From Ref. [179].

Assuming an isotropic tau neutrino flux, the number of taus that emerge from the Earth with
sufficient energy to be detected is proportional to an "effective solid angle"

fi e ff = J P{9) cos 9 d cos 9 d(\.
Evaluation of the integrals [286] yields the unfortunate expression

fLff — 27T-

T v

3^ T /^ot _ 1

2R(n

(201)

(202)

At the relevant energies, however, the neutrino interaction length satisfies L\ ot <C R®. In addition, if the
hypothesized non-perturbative cross section enhancement is less than typical hadronic cross sections, we
have L" ot » L T . With these approximations, Eq. (202) simplifies to [287]

^eff

tv2tt

2n LtotL

4i?|L£ c

(203)

Equation (203) describes the functional dependence of the Earth- skimming event rate on the non-
perturbative cross section. This rate is, of course, also proportional to the neutrino flux $^ a11 at Eq. Thus,
the number of Earth- skimming neutrinos is given by

Nes ~ C E s

Cj^all

T CC

(204)

{(T cc -h cr Np j

where Ces is the number of Earth- skimming events expected for some benchmark flux <3>Q a11 in the
absence of new physics.

In contrast to Eq. (204), the rate for quasi-horizontal showers has the form

Nqh = Cqh

(205)

74

Fig. 35: Event rates for Earth- skimming (left) and quasi-horizontal (right) events in the $^ all /$Q a11 — ctnp/o"cc
plane. Note that the contours are roughly orthogonal, and so the two types of event provide complementary
information about flux and cross section.

Fig. 36: Projected determination of neutrino fluxes and cross sections at y/s ^ 250 TeV from future Auger data.
The different shaded regions indicate the 90%, 95%, 99% and 3a confidence level contours in the $^ all /$Q a11 —
ctnp/ctcc plane, for 7V°jf = 1, N$g = 10 (left), N°f = 1, N$g = 7 (middle), and N°f = 1, = 5 (right).
The dashed line indicates the result of including the systematic uncertainty on the NLO QCD CC neutrino-nucleon
cross section. FromRef. [294].

where Cqh is the number of quasi-horizontal events expected for flux $Q a11 .

Now, without loss of generality, we normalize the neutrino flux to the Waxmann-Bahcall bound,
i.e. <3>Q a11 = ^wb- We use ^ e ex P ecte d rates f° r the the benchmark flux to determine the values of
Ces and Cqh in Eqs. (204) and (205), (Ces = 0.15 and Cqh = 0.06, as shown in Table 7). Given a
flux $^ a11 and new non-perturbative physics cross section a^ p , both TVes and A/qh are determined. On
the other hand, given just a quasi-horizontal event rate A/qh, it is impossible to differentiate between an
enhancement of the cross section due to non-perturbative physics and an increase of the flux. However, in
the region where significant event rates are expected, the contours of A/qh and N^s, given by Eqs. (204)
and (205), are more or less orthogonal and provide complementary information. This is illustrated in
Fig. 35. With measurements of Nq^ and N^ s , both and $^ a11 may be determined independently,
and neutrino interactions beyond the (perturbative) SM may be unambiguously identified.

75

We now turn to determining the projected sensitivity of Auger to neutrino fluxes and cross sections.
The quantities TVes and TVqh as defined in Eqs. (204) and (205) can be regarded as the theoretical values
of these events, corresponding to different points in the cfr^aii/cfr^aii _ a^p/vcc parameter space. For a
given set of observed rates N^ s and Nq^, two curves are obtained in the two-dimensional parameter
space by setting A^jf = A^es and Nq^ = Aqh- These curves intersect at a point, yielding the most
probable values of flux and cross section for the given observations. Fluctuations about this point define
contours of constant x 2 in an approximation to a multi-Poisson likelihood analysis [304]. The contours
are defined by

where i = ES, QH [305]. In Fig. 36, we show results for three representative cases. Assuming (iVgg s = 1,
Aqjj = 10), (N°f = 1, iVgg = 7)? and (7V obs = 1? jyobs = 5) we show the 90%9 95%9 99% and 3a C L

contours for 2 d.o.f. (x 2 = 4.61, 5.99, 9.21, and 11.83, respectively). For N°f = 1 and A^ s = 10, the
possibility of a SM interpretation along the a^ p = axis (taking into account systematic uncertainties)
would be excluded at greater than 99% CL for any assumed flux.

In summary, we have found that observation of 1 Earth- skimming and 10 quasi-horizontal events
would exclude the (perturbative) SM at the 99% CL. Thus the expected low neutrino "luminosity" is not
at all a show-stopper, and the observatory has the potential to uncover new physics at scales exceeding
those accessible to the LHC. If new non-perturbative physics exists, a decade or so would be required to
uncover it at Auger in the best case scenario (cosmic neutrino flux at the Waxman-Bahcall level and vN
cross section about an order of magnitude above the SM prediction). Of course future CR experiments
should benefit from a much larger aperture, making such discovery conceivable in a short time scale. An
optimist may even imagine the possibility of probing QCD dynamics [306].

Acknowledgments

I would like to thank Felix Aharonian, Markus Ahlers, Segev BenZvi, Hans Bliimer, Mandy Cooper
Sarkar, Olivier Deligny, Tere Dova, Jonathan Feng, Haim Goldberg, Concha Gonzalez Garcia, Yann
Guardincerri, Francis Halzen, Nicole Krohm, Fred Kuehn, Antoine Letessier Selvon, Jim Matthews,
Teresa Montaruli, Giulia Pancheri, Tom Paul, Subir Sarkar, Viviana Scherini, Al Shapere, Paul Sommers,
John Swain, Diego Torres, Alan Watson, and Tom Weiler for valuable discussions and permission to
reproduce some of the figures. I am grateful to Stefanie Pinnow for her proofreading skills and to my
sister Gaby for her assistance with slide management. L.A.A. is supported by the U.S. National Science
Foundation (NSF Grant No PHY-0757598) and the UWM Research Growth Initiative. Any opinions,
findings, and conclusions or recommendations expressed in this material are those of the author and do
not necessarily reflect the views of the National Science Foundation.

Appendices

A Cosmogenic /3-DK and A* processes

If UHECRs are heavy nuclei, the relic photons can excite the GDR at nucleus energies E > 10 11 GeV
(corresponding to 10 MeV — 30 MeV in the nuclear rest frame), and thus there should be

accompanying photo-dissociated free neutrons, themselves a source of /3-decay antineutrinos. The decay
mean free path of a neutron is cY n r n = (E n /10 u GeV) Mpc, the lifetime being boosted from its rest
frame value r n = 886 s to its lab value via T n = E n /m n . Compared to cosmic distances > 100 Mpc,
the decay of even the boosted neutron may be taken as nearly instantaneous, and thus all free neutrons are
themselves a source of /3-decay cosmogenic antineutrinos. The neutron emissivity Q n (E n ), defined as
the mean number of neutrons emitted per comoving volume per unit time per unit energy as measured at
the source can be estimated as follows. Neutrons with energies above 10 9 3 GeV have parent iron nuclei
with r > Tq ~ 2 x 10 9 which are almost completely disintegrated in distances of less than 100 Mpc (see

(206)

76

Sec. 2.3.2). Thus, it is reasonable to define a characteristic time r r given by the moment at which the
number of nucleons is reduced to 1 je of its initial value A, and presume the nucleus, emitted at distance d
from the Earth, is a traveling source that at D « (d — cr r ) disintegrates into A nucleons all at once [120].
Then, the number of neutrons with energy E n = Ea/A can be approximated by the product of 1/2 the
number of nucleons generated per nucleus and the number of nuclei emitted, i.e. Q n (E n ) =NQa, where
N = A — Z is the mean neutron number of the source nucleus. Now, to obtain an estimate of the
diffuse antineutrino flux one needs to integrate over the population of nucleus-emitting-sources out to
the horizon [307, 308]

1 f 1T , ^ ^ ^ \ ( Dm n \] f Q 7 dP , N

tio J |_ \ E n r n J ] J dev

6 [E^- E n e 17 (l + cos O^/rrin] , (A.l)

47T

^ dcosOjj

where the r 2 in the volume integral is compensated by the usual 1/r 2 fall-off of flux per source. Here,
Hq is the Hubble constant, E v and E n are the antineutrino and neutron energies in the lab, 6 V is the
antineutrino angle with respect to the direction of the neutron momentum in the neutron rest-frame,
and e„ is the antineutrino energy in the neutron rest-frame. The last three variables are not observed
by a laboratory neutrino-detector, and so are integrated over. The observable E v is held fixed. The
delta- function relates the neutrino energy in the lab to the three integration variables. The parameters
appearing in Eq. (A.l) are the neutron mass and rest-frame lifetime (m n and r n ). Finally, dP/dejj
is the normalized probability that the decaying neutron produces a V with energy in the neutron
rest-frame. Note that the maximum V energy in the neutron rest frame is very nearly Q = m n —
m p — m e = 0.71 MeV. Integration of Eq. (A.l) can be easily accomplished, especially when two good
approximations are applied [309] . The first approximation is to think of the /3-decay as a 1 —> 2
process of 5m n -> e~ + 17, in which the neutrino is produced monoenergetically in the rest frame, with
ej7 = eo ~ SrriN(l — m 2 e / 5 2 rriN) /2 ~ 0.55 MeV, where Stun ~ 1.30 MeV is the neutron-proton mass
difference. In the lab, the ratio of the maximum V energy to the neutron energy is 2eo/m n ~ 10 -3 ,
and so the boosted V has a spectrum with E v e (0, 10 -3 E n ). The second approximation is to replace
the neutron decay probability 1 - e - Dm n/E n r n w ^ a ste p f unc ti n 6(£^f ax - E n ) at some energy
^max ^ q^jj mn /r n ) = (D/10 Mpc) x 10 12 GeV Combining these two approximations we obtain

where E™ m = max{^ ; r , A ^ Ev }, and £^™ ax is the energy cutoff at the nucleus-emitting-source
< A(D/W Mpc) x 10 12 GeV. For Q A oc E^ a , integration of Eq. (A.2) leads to

Ea,f

^ e (Eu) w 10 6
where Ep > 10 6 3 (56/A) GeV and

ZTTnax
A

rmax \ csl

The sub-PeV antineutrino spectrum is flat, as all the free neutrons have sufficient energy E n > Er /A,
to contribute equally to all the V energy bins below 10 6 GeV. Taking a = 2 as a reasonable example,

and inputting the observational value for iron nuclei, E\ ro Jqr

becomes

w 10 5 eV m~ 2 s" 1 sr" 1 , Eq. (A.3)

ro

El <& Ve w 4 x 10 1 ( ^ ) eVm" 2 s" 1 sr 1 . (A.5)

77

Note that the /3-decay process gives initial antineutrino flavor ratios 1:0:0 and Earthly ratios nearly 3 :
1:1. Compared to full-blown Monte Carlo simulations [310,311], this back of the envelope calculation
underestimates the flux by about 30%. This is because the preceeding calculation does not account for
possible neutrons produced in £>7cmb collisions. Of course the situation described above represents the
most extreme case, in which all cosmic rays at the end of the spectrum are heavy nuclei. A more realistic
guess would be that the composition at the end of the spectrum is mixed.

Photodisintegration of high energy nuclei is followed by immediate photo-emission from the ex-
cited daughter nuclei. For brevity, we label the photonuclear process A + 7 — > A'* + X, followed by
A f * —> A' + 7-ray as "A*" [312]. The GDR decays dominantly by the statistical emission of a single
nucleon, leaving an excited daughter nucleus (A — 1)*. The excited daughter nuclei typically de-excite
by emitting one or more photons of energies e^ xn ^1 — 5 MeV in the nuclear rest frame. The lab-frame
energy of the 7-ray is then ei^ AB = Ta ^ xn , where Ta = Ea/ttia is the boost factor of the nucleus in
the lab.

Of interest here is the 7-ray flux produced when the photo-dissociated nuclear fragments produced
on the CMB and CIB de-excite. These 7-rays create chains of electromagnetic cascades on the CMB and
CIB, resulting in a transfer of the initial energy into the so-called Fermi-LAT region, which is bounded
by observation to not exceed cj cas ~ 5.8 x 10 _7 eV/cm 3 [170]. Fortunately, we can finesse the details
of the calculation by arguing in analogy to work already done. The photodisintegration chain produces
one /3-decay neutrino with energy of order 0.5 MeV in the nuclear rest frame, for each neutron produced.
Multiplying this result by 2 to include photodisintegration to protons in addition to neutrons correctly
weights the number of steps in the chain. Each step produces on average one photon with energy ~
3 MeV in the nuclear rest frame. In comparison, about 12 times more energy is deposited into photons.
Including the factor of 12 relating uo 1 to ojp e , we find from (A.5) that cosmogenic photodisintegration/de-
excitation energy, cj cas ^ 1.4 x 1CT 11 (56/A) eV/cm 3 , is more than three orders of magnitude below
the Fermi-LAT bound. 10 This result appears to be nearly invariant with respect to varying the maximum
energy of the Fe injection spectrum (with a larger £ max , the additional energy goes into cosmogenic pion
production). Thus, there is no constraint on a heavy nuclei cosmic-ray flux from the A* mechanism.

Neutron emission from a cosmological distribution of CR sources would also lead to a flux of V e .
In analogy to our previous estimate, we sum over the neutron-emitting sources out to the edge of the
universe; Eq. (A.2) becomes [314]

where Q n {E n ) is the neutron volume emissivity. An upper limit can be placed on Q n by assum-
ing an extreme situation in which all the CRs escaping the source are neutrons, i.e., ej^ > 10 1 =
J dE n E n Q n (E n ). With the production rate of ultrahigh energy protons (91), and an assumed injection
spectrum Q n oc E~ 2 , Eq. (A.6) gives

which is about three orders of magnitude below the Waxman-Bahcall bound.
B Energy density of the EM cascade

EM interactions of photons and leptons with the extragalactic radiation backgrounds and magnetic field
can happen on time-scales much shorter than their production rates. The relevant processes with back-
ground photons contributing to the differential interaction rates 7 ee , 7 7e , and 7 e7 are inverse Comp-
ton scattering (ICS), e ± + 7b gr —> e ± + 7, pair production (PP), 7 + 7b gr —> e + + e~, double pair

10 These rough calculations are consistent with the more rigurous analysis presented in [313].

(A.6)

El ^ w 3 x 10" 11 GeVcm^s" 1 sr

-1

(A.7)

78

production (DPP) 7 + 7b gr —> e + + e + e + + e , and triple pair production (TPP), e ± + 7b gr —>
e ± + e+ + e~ [106,315].

High energetic electrons and positrons may also lose energy via synchrotron radiation in the in-
tergalactic magnetic field B with a random orientation sin 9 with respect to the velocity vector. We will
assume in the following that the field strength B = 10 -12 G. 11 This leads to an efficient transfer of
energy into the EM cascade. The synchrotron power spectrum (W eV _1 ) has the form

P(£ e , £? 7 ) = F(Ej/E c ) ; F(t) = i / <fe tf 5/3 (z) ,

Z7T m e j t

(B.l)

where K 5 / 3 is the modified Bessel function and E c = (3e J Bsin0/2m e )(^ e /m e ) 2 . This can be treated
as a continuous energy loss of the electrons and positrons with a parameter 12

W*e) = ^ / decs* / ^ 7 W,£ 7 ) = f (^) (^) • (B.2)

We will assume in the calculation that the intergalactic magnetic field is primordial with a (flux-conserving)
redshift dependence B(z) = (1 + z) 2 B(0). Note that the synchrotron energy loss has then a redshift
dependence similar to BH pair production in the CMB, i.e. b syn (z, E) = (1 + ^) 2 6 syn (0, (1 + z)E). It is
also convenient to define 7^ n (i? e , E 7 ) = V(E ei E 1 )/E 1 , which has an analogous redshift dependence,
i.e. ffi(z, E e , = (1 + z) 4 7 !7 n (0> (1 + z)E e , (1 + z)E 1 ).

The fast evolution of the cascade is governed by the set of differential equations,

d i Y e {E)=-T e (E)Y e {E) + d E {b{E)Y e {E)) + J dE' [ lie {E' , E)Y 1 {E') + lee (E' , E)Y e {E')] ,

d t Yy(E) = - T 7 (E)Y 7 (E) + J dE , V ^ E ^ Y e (E') + J dE'^ ei (E r , E)Y e (E r ) , (B.3)

which determines the evolution on short time-scales At Y pi <C 1 (the redshift z is kept fixed meanwhile).
The initial condition Yy/ e (E)\f =0 is given by the sum of previously developed cascades and the newly
generated contributions from proton interactions.

The solution of the system (B.3) for an infinitesimally small step At can be written for a discrete
energy spectrum, Ni ~ AEiYi, as

With the transition matrix T(At), defined by Eq. (B.4), we can efficiently follow the development of the
EM cascade over a distance At = 2 N At via matrix doubling [167]:

T{2 N At)~ [T(At)] 2N . (B.5)
The total energy of the cascade can be obtained by integrating Eq. (120):

J dE E^ s (z, E) = - J dEEdg [6 cas (z, S)2t v {z, E)\ . (B.6)

Integrating the r.h.s. by parts yields

r.h.s. = - J dEd E E^-b^{z,S)^ v {z,E) + J dE^-b cas (z^)^ p (z, E) . (B.7)

n For B > 10~ 12 G, TPP by electrons can be neglected [316]. For E < 10 12 GeV, DPP of photons can also be safely
neglect in the calculation [317].

12 Note the identity f dE [E fe(6n e ) + f dE'V(E\ E)n e ] = 0, implying overall energy conservation.

d
dt

79

io- 2 0.1 1 10 10 2 10 3 10 4 10 5 10 6 10 7 10 8 10 9 10 10 10 n 10 12
E [GeV]

Fig. B.l: Cosmogenic photon (blue lines) and neutrino (red line) fluxes from a "vanilla" proton test- spectrum
(shown as black line) of UHECRs, assuming two exponential cutoffs (E m i n = 10 8 GeV and E max = 10 11,5 GeV),
power index 7 = 2 and evolution with n = 3 and < z < 2. The fluxes are computed for 3 extragalactic magnetic
field strengths: 1 pG, 1 nG, and 50 nG. Proton diffusion in extragalactic magnetic fields B » nG has not yet been
included in the calculation. We have neglected radiative losses of muons and pions before decaying into neutrinos;
the energy loss in muons is down by a factor 10 6 compared to electrons. This figure is courtesy of Markus Ahlers.

The first term vanishes since 6 cas = for sufficiently low energies and 3? v = beyond the maximal
energy. The time integration of the l.h.s. between the present epoch (t = 0) and the first sources (t max )
gives

J dt[l.h.s.] = J dEEn cas (E) = u cas , (B.8)

hence we obtain Eq. (121).

Synchrotron radiation in strong magnetic fields, 1 nG < B < 50 nG, can also by-pass the EM cas-
cade and transfer energy to sub-GeV photons that are unconstrained by the Fermi-LAT spectrum [318].
This would be relevant for electrons around E e ~ 10 9 GeV, where synchrotron loss starts to dominate
over ICS loss in the CMB; see Fig. B.l. From Eqs. (20) and (27) we can verify that (E* yn ) < 1 GeV. Of
course, reaching a 50 nG field requires some astrophysical generator to augment the primordial field. 13
For simplicity, we have taken the redshift evolution of the generator to be the same as the evolution of
the primordial field. If were to assume an evolution which follows SFR, this will not significantly alter
our conclusions.

C TOTal Elastic and diffractive cross section Measurement

The TOTEM experiment is dedicated to the measurement of the total pp cross section with the luminosity
independent method based on the Optical Theorem,

a tot = Im/(0) , (C.l)

where /(#) satisfies

at s ail s

13 For Xb ~ 1 Mpc, the current upper limit on primordial seed fields is about 1 nG [319, 320].

80

with i) the angle of the scattered proton with respect to the beam direction. Squaring Eq. (C.l) we obtain

2 16tt Im 2 /(0) da el

a i —

'tot

Re 2 /(0) + Im 2 /(0) dt

16tt [dN el /dt] t=0 _

t=o

1 + P 2 C

where C is the integrated luminosity. Now, following [321, 322], we can obtain the total cross section
independently from C, by using a tot = (a e \ + cr ine i) = (N e \ + 7V ine i)/£,

167T [dN el /dt] t=0

1 + p 2 (N el + 7V ine i)

Here, N e \ and N[ ne \ are the numbers of elastic and inelastic events, and p = 0.10 ± 0.01 is the ratio
between the real and imaginary parts of the forward scattering amplitude [323]. 14 The difficult aspect
of this measurement is obtaining a good extrapolation of the cross section for low momentum transfer.
Recall that —t = s(l — cos #)/2 ~ si} 2 / 4, and so t — >> implies a measurement in the extreme forward
direction. The TOTEM experiment aims to measure down to values of \t\ ~ xlO -3 GeV 2 , which
corresponds to i? « 4.5 /irad [324].

D Raiders of the lost holy grail

In 1976 't Hooft observed that the Standard Model does not strictly conserve baryon and lepton num-
ber [325, 326]. Rather, non-trivial fluctuations in SU{2) gauge fields generate an energy barrier interpo-
lating between topologically distinct vacua. An index theorem describing the fermion level crossings in
the presence of these fluctuations reveals that neither baryon nor lepton number is conserved during the
transition, but only the combination B — L. Inclusion of the Higgs field in the calculation modifies the
original instanton configuration [327]. An important aspect of this modification (called the "sphaleron")
is that it provides an explicit energy scale of about 10 TeV for the height of the barrier. This barrier can
be overcome through thermal transitions at high temperatures [328-330], providing an important input to
any calculation of cosmological baryogenesis. More speculatively, it has been suggested [331-333] that
the topological transition could take place in two particle collisions at very high energy. The anomalous
electroweak contribution to the partonic process can be written as

di{s) = 5.3 x 10 3 mb • e ~^ /aw) Fw ^ e) , (D.l)

where aw — 1/30, the tunneling suppression exponent i^V(e) is sometimes called the "holy-grail func-
tion", and e = >/§ / \\.i\mw I otw) — V^/30 TeV. Thus, it is even possible that at or above the sphaleron
energy the cross-section could be of (9(mb) [334]. Of particular interest here would be an enhancement
of the neutrino cross section over the perturbative SM estimates, say by an order of magnitude in the
energy range 9.5 < log 10 (^/GeV) < 10.5.

The argument for strong damping of the anomalous cross section for y/§ > 30 TeV was convinc-
ingly demonstrated in [335, 336], in the case that the classical field providing the saddle point interpo-
lation between initial and final scattering states is dominated by spherically symmetric configurations.
This 0(3) symmetry allows the non- vacuum boundary conditions to be fully included in extremizing
the effective action. In [337] it was shown that a sufficient condition for the 0(3) dominance is that the
interpolating field takes the form of a chain of "lumps" which are well- separated, so that the each lump
lies well into the exponentially damped region of its nearest neighbors. However, we are not aware of
any reason that such lumped interpolating fields should dominate the effective action. It is thus of interest
to explore the other extreme, in which non-spherically symmetric contributions dominate the effective
action (and let experiment be the arbiter).

14 Note that the quoted value of p is an extrapolation to = 14 TeV, and may be measured by the LHC experiments.
Otherwise, it will contribute to the uncertainty in cr to t.

81

It was shown [334] that for the simple sphaleron configuration s-wave unitarity is violated for
y/§ > 4ttMw/(^w ~ 36 TeV. If for y/§ > 36 TeV we saturate unitarity in each partial wave, then this
yields a geometric parton cross section ttR 2 , where R is some average size of the classical configuration.
As a fiducial value we take the core size of the Manton-Klinkhamer sphaleron, R ~ 4/Mw — 10 -15 cm.
In this simplistic model, the uN cross section is

where x min = s min /s = (36) 2 /2ME u ~ 0.065. In the region 0.065 < x < 3(0.065) the PDF for the
up and down quarks is well approximated by / ~ 0.5 /x, so the expression for the cross section becomes

where the last factor of 2/2 takes into account the (mostly) 2 contributing quarks (u, d) in this range of
x, and the condition that only the left-handed ones contribute to the scattering. This is about 80 times
the SM cross section. Of course this calculation is very approximate and the cross section can easily be
smaller by a factor of 10 {e.g. if R is 1/3 of the fiducial value used).

References

[1] V. F. Hess, Phys. Z. 13, 1804 (1912).

[2] P. Auger, R. Maze, T. Grivet-Meyer, Comptes Rendus 206, 1721 (1938).

[3] P. Auger, P. Ehrenfest, R. Maze, J. Daudin, Robley, and A. Freon, Rev. Mod. Phys. 11, 288 (1939).
[4] M. Nagano and A. A. Watson, Rev. Mod. Phys. 72, 689 (2000).
[5] A. Letessier-Selvon and T. Stanev, arXiv: 1103.0031 [astro-ph.HE].

[6] L. Anchordoqui, T. C. Paul, S. Reucroft and J. Swain, Int. J. Mod. Phys. A 18, 2229 (2003)

[arXiv:hep-ph/0206072].
[7] T. Abu-Zayyad et al. 9 Nucl. Instrum. Meth. A 450, 253 (2000).
[8] R. M. Baltrusaitis et al., Nucl. Instrum. Meth. A 240, 410 (1985).
[9] J. Abraham et al. [Pierre Auger Collaboration], Nucl. Instrum. Meth. A 523, 50 (2004).
[10] J. Abraham et al. [Pierre Auger Collaboration], Nucl. Instrum. Meth. A 613, 29 (2010).
[11] J. A. Abraham et al. [Pierre Auger Collaboration], Nucl. Instrum. Meth. A 620, 227 (2010)
[arXiv:0907.4282].

[12] P. Abreu et al. [Pierre Auger Collaboration], Astropart. Phys. 34, 368 (2011) [arXiv: 1010.6162
[astro-ph.HE]].

[13] Y. Takahashi et al. [JEM-EUSO Collaboration], New J. Phys. 11, 065009 (2009) [arXiv:0910.4187

[astro-ph.HE]].
[14] T. Weiler, private communication.

[15] C. Amsler et al. [Particle Data Group], Phys. Lett. B 667, 1 (2008).
[16] M. Aglietta et al. [EAS-TOP Collaboration], Astropart. Phys. 21, 583 (2004).
[17] T. Antoni et al. [KASCADE Collaboration], Astropart. Phys. 24, 1 (2005) [arXiv:astro-
ph/0505413].

[18] W. D. Apel et al [KASCADE Collaboration], Astropart. Phys. 31, 86 (2009) [arXiv:08 12.0322
[astro-ph]].

[19] J. Bluemer, R. Engel and J. R. Hoerandel, Prog. Part. Nucl. Phys. 63, 293 (2009) [arXiv:0904.0725
[astro-ph.HE]].

!min partons

(D.2)

<lackdisk(^) - 7ri? 2 (0.5) (In 3) (2/2)

~ 1.5 x 10" 30 cm 2

(D.3)

82

[20] M. Hillas, Phys. Lett. A 24, 677 (1967).

[21] V. Berezinsky, A. Z. Gazizov and S. I. Grigorieva, Phys. Rev. D 74, 043005 (2006) [arXiv:hep-
ph/0204357].

[22] R. Aloisio, V. Berezinsky, P. Blasi, A. Gazizov, S. Grigorieva and B. Hnatyk, Astropart. Phys. 27,

76 (2007) [arXiv:astro-ph/0608219].
[23] G. Cocconi, Astropart. Phys. 4, 281 (1996).
[24] A. A. Penzias and R. W. Wilson, Astrophys. J. 142, 419 (1965).
[25] K. Greisen, Phys. Rev. Lett. 16, 748 (1966).

[26] G. T. Zatsepin and V. A. Kuzmin, JETP Lett. 4, 78 (1966) [Pisma Zh. Eksp. Teor. Fiz. 4, 114
(1966)].

[27] J. N. Bahcall and E. Waxman, Phys. Lett. B 556, 1 (2003) [arXiv:hep-ph/0206217].
[28] R. U. Abbasi et al. [HiRes Collaboration], Phys. Rev. Lett. 100, 101101 (2008) [arXiv:astro-
ph/0703099].

[29] J. Abraham et al. [Pierre Auger Collaboration], Phys. Rev. Lett. 101, 061101 (2008)
[arXiv:0806.4302].

[30] J. Abraham et al. [Pierre Auger Collaboration], Phys. Lett. B 685, 239 (2010) [arXiv: 1002.1975
[astro-ph.HE]].

[31] J. Linsley and A. A. Watson, Phys. Rev. Lett. 46, 459 (1981).

[32] L. D. Landau and I. Pomeranchuk, Dokl. Akad. Nauk Ser. Fiz. 92, 535 (1953).

[33] A. B. Migdal, Phys. Rev. 103, 1811 (1956).

[34] J. Abraham et al. [Pierre Auger Collaboration], Astropart. Phys. 27, 155 (2007) [astro-ph/0606619].
[35] J. Abraham et al. [Pierre Auger Collaboration], Astropart. Phys. 29, 243 (2008) [arXiv:0712.1 147
[astro-ph]].

[36] J. Abraham et al. [Pierre Auger Collaboration], Astropart. Phys. 31, 399 (2009) [arXiv:0903.1 127
[astro-ph.HE]].

[37] R. U. Abbasi et al. [HiRes Collaboration], Phys. Rev. Lett. 104, 161101 (2010) [arXiv:09 10.4 184
[astro-ph.HE]].

[38] D. R. Bergman [HiRes Collaboration], Nucl. Phys. Proc. Suppl. 136, 40 (2004) [arXiv:astro-
ph/0407244].

[39] D. J. Bird et al. [HiRes Collaboration], Phys. Rev. Lett. 71, 3401 (1993).

[40] J. Abraham et al. [Pierre Auger Collaboration], Phys. Rev. Lett. 104, 091101 (2010)

[arXiv: 1002.0699 [astro-ph.HE]].
[41] R. Ulrich, R. Engel, S. Muller, F. Schussler and M. Unger, Nucl. Phys. Proc. Suppl. 196, 335 (2009)

[arXiv:0906.3075 [astro-ph.HE]].
[42] L. A. Anchordoqui, M. T. Dova, L. N. Epele and S. J. Sciutto, Phys. Rev. D 59, 094003 (1999)

[arXiv:hep-ph/9810384].
[43] J. Linsley, Phys. Rev. Lett. 34, 1530 (1975).

[44] P. Sommers, Astropart. Phys. 14, 271 (2001) [arXiv:astro-ph/0004016].
[45] P. Sokolsky, P. Sommers and B. R. Dawson, Phys. Rept. 217, 225 (1992).
[46] J. Aublin, E. Parizot, Astron. Astrophys. 441, 407 (2005).
[47] R. Bonino et al. , in preparation.

[48] P. Abreu et al. [Pierre Auger Collaboration], Astropart. Phys. 34, 267 (2011). [arXiv: 1103.2721
[astro-ph.HE]].

[49] N. Hayashida et al. [AGASA Collaboration], Astropart. Phys. 10, 303 (1999) [arXiv:astro-
ph/9807045].

83

[50] T. Stanev, Astrophys. J. 479, 290 (1997) [arXiv:astro-ph/9607086].
[51] J. Candia, S. Mollerach and E. Roulet, JCAP 0305, 003 (2003) [arXiv:astro-ph/0302082].
[52] A. Calvez, A. Kusenko and S. Nagataki, Phys. Rev. Lett. 105, 091101 (2010) [arXiv: 1004.2535
[astro-ph.HE]].

[53] A. H. Compton and I. A. Getting, Phys. Rev. 47, 817 (1935).

[54] M. Kachelriess and P. D. Serpico, Phys. Lett. B 640, 225 (2006) [arXiv:astro-ph/0605462].
[55] J. Abraham et al. [Pierre Auger Collaboration], Science 318, 938 (2007) [arXiv:07 1 1 .2256 [astro-
ph]].

[56] J. Abraham et al. [Pierre Auger Collaboration], Astropart. Phys. 29, 188 (2008) [Erratum-ibid. 30,

45 (2008)] [arXiv:07 12.2843 [astro-ph]].
[57] R. U. Abbasi etal, Astropart. Phys. 30, 175 (2008) [arXiv:0804.0382 [astro-ph]].
[58] P. Abreu et al. [Pierre Auger Collaboration], Astropart. Phys. 34, 314 (2010) [arXiv: 1009. 1855].
[59] D. S. Gorbunov, P. G. Tinyakov, 1. 1. Tkachev, S. V. Troitsky, [arXiv:0804.1088 [astro-ph]].
[60] A. M. Hillas, Ann. Rev. Astron. Astrophys. 22, 425 (1984).
[61] W. F. G. Swann, Phys. Rev. 43, 217 (1933).
[62] E. Fermi, Phys. Rev. 75, 1169 (1949).

[63] T. K. Gaisser, Cosmic Rays And Particle Physics, (Cambridge, UK: Univ. Press, 1990).
[64] R. Blandford and D. Eichler, Phys. Rept. 154, 1 (1987).

[65] M. Ahlers, L. A. Anchordoqui, J. K. Becker, T. K. Gaisser, F. Halzen, D. Hooper, S. R. Klein. P.

Meszaros, S. Razzaque, and S. Sarkar, FERMILAB-FN-0847-A, YITP-SB- 10-01.
[66] L. O. Drury, Contemp. Phys. 35, 231 (1994).
[67] A. R. Bell, Mon. Not. Roy. Astron. Soc. 182, 147 (1978).
[68] R. D. Blandford and J. P. Ostriker, Astrophys. J. 221, L29 (1978).

[69] L. A. Anchordoqui, G. E. Romero and J. A. Combi, Phys. Rev. D 60, 103001 (1999) [arXiv:astro-
ph/9903145].

[70] D. F. Torres and L. A. Anchordoqui, Rept. Prog. Phys. 67, 1663 (2004) [arXiv:astro-ph/0402371].
[71] B. L. Fannaroff and J. M. Riley, Mon. Not. Roy. Astron. Soc. 167, 31 (1974).
[72] R. D. Blandford and M. J. Rees, Mon. Not. Roy. Astron. Soc. 169, 395 (1974).
[73] A. Rosen, P. A. Hughes, G. C. Duncan and P. E. Hardee, Astrophys. J. 516, 729 (1999) [arXiv:astro-
ph/9901046].

[74] M. C. Begelman, R. D. Blandford and M. J. Rees, Rev. Mod. Phys. 56, 255 (1984).

[75] L. O. Drury, Rept. Prog. Phys. 46, 973 (1983).

[76] A. N. Kolmogorov, C. R. Acad. URSS 30, 201 (1941).

[77] J. P. Rachen and P. L. Biermann, Astron. Astrophys. 272, 161 (1993) [arXiv:astro-ph/9301010].
[78] G. B. Rybicki and A. P. Lightman, Radiative Processes in Astrophysics (New York: Wiley-

Interscience, 1979).
[79] T. A. Armstrong et al, Phys. Rev. D 5, 1640 (1972).
[80] P. L. Biermann and P. A. Strittmatter, Astrophys. J. 322, 643 (1987).

[81] L. A. Anchordoqui and H. Goldberg, Phys. Rev. D 65, 021302 (2002) [arXiv:hep-ph/0106217].
[82] F. P. Israel, Astron. Astrophys. Rev. 8, 237 (1998).

[83] J. J. Condon, G. Helou, D. B. Sanders, and B. T. Soifer, Astrophys. J. Suppl. 103, 81 (1996).
[84] F. Aharonian [H.E.S.S. Collaboration], Astrophys. J. 695, L40 (2009) [arXiv:0903.1582 [astro-
ph.CO]].

[85] M. J. Hardcastle, D. M. Worrall, R. P. Kraft, W. R. Forman, C. Jones and S. S. Murray, Astrophys.
J. 593, 169 (2003) [arXiv:astro-ph/0304443].

84

[86] J. O. Burns, E. D. Feigelson, and E. J. Schreier, Astrophys. J. 273, 128 (1983).

[87] G. E. Romero, J. A. Combi, L. A. Anchordoqui and S. E. Perez Bergliaffa, Astropart. Phys. 5, 279
(1996) [arXiv:gr-qc/9511031].

[88] M. Honda, Astrophys. J. 706, 1517 (2009) [arXiv:09 11.0921 [astro-ph.HE]].

[89] N. Junkes, R. F. Haynes, J. I. Harnett, and D. L. Jauncey, Astron. Astrophys. 269, 29 (1993) [Erra-
tum, ibid 274, 1009(1993)].

[90] P. Sreekumar, D. L. Bertsch, R. C. Hartman, P. L. Nolan and D. J. Thompson, Astropart. Phys. 11,
221 (1999) [arXiv:astro-ph/9901277].

[91] J. E. Grindlay et al, Astrophys. J. 197, L9 (1975).

[92] A. A. Abdo et al. [Fermi-LAT Collaboration], Science 328, 725 (2010) [arXiv: 1006.3986 [astro-
ph.HE]].

[93] V. L. Ginzburg and S. I. Syrovatskii, Ann. Rev. Astron. Astrophys. 3, 297 (1965).

[94] G. J. Fishman and C. A. Meegan, An. Rev. Astron. Astrophys. 33, 415 (1995).

[95] B. Link and R. I. Epstein, Astrophys. J. 466, 764 (1996) [arXiv:astro-ph/9601033].

[96] C. A. Meegan et. al., Nature 355, 143 (1992).

[97] M. R. Metzger et. al., Nature 387, 878 (1997).

[98] T. Piran, Phys. Rept. 314, 575 (1999) [arXiv:astro-ph/9810256].

[99] P. Meszaros, Rept. Prog. Phys. 69, 2259 (2006) [arXiv:astro-ph/0605208].

[100] E. Waxman, Lect. Notes Phys. 598, 393 (2003). [arXiv:astro-ph/0303517].

[101] E. Waxman, Phys. Rev. Lett. 75, 386 (1995) [arXiv:astro-ph/9505082].

[102] W. Coburn and S. E. Boggs, Nature 423, 415 (2003).

[103] E. Nakar, T. Piran and E. Waxman, JCAP 0310, 005 (2003) [arXiv:astro-ph/0307290].
[104] F. W. Stecker, Phys. Rev. Lett. 21, 1016 (1968).

[105] V. S. Berezinsky and S. I. Grigor'eva, Astron. Astrophys. 199, 1 (1988).

[106] G. R. Blumenthal, Phys. Rev. D 1, 1596 (1970).

[107] R. M. Barnett et al. [Particle Data Group], Phys. Rev. D 54, 1 (1996).

[108] L. Montanet et al. [Particle Data Group], Phys. Rev. D 50, 1 173 (1994). See p. 1335.

[109] I. Golyak, Mod. Phys. Lett. A 7, 2401 (1992).

[110] L. A. Anchordoqui, M. T. Dova, L. N. Epele and J. D. Swain, Phys. Rev. D 55, 7356 (1997)
[arXiv:hep-ph/9704387].

[Ill] M. Abramowitz, I. A. Stegun, "Handbook of Mathematical Functions" (Dover, New York, 1970).

[112] J. L. Puget, F W. Stecker and J. H. Bredekamp, Astrophys. J. 205, 638 (1976).

[113] M. J. Chodorowski, A. A. Zdziarski, and M. Sikora, Astrophys. J. 400, 181 (1992).

[114] S. Michalowski, D. Andrews, J. Eickmeyer, T. Gentile, N. Mistry, R. Talman and K. Ueno, Phys.

Rev. Lett. 39, 737 (1977).
[115] F. W. Stecker, Phys. Rev. 180, 1264 (1969).
[116] E. Hayward, Rev. Mod. Phys. 35, 324 (1963).

[117] M. H. Salamon and F. W. Stecker, Astrophys. J. 493, 547 (1998) [arXiv:astro-ph/9704166].

[118] L. N. Epele and E. Roulet, JHEP 9810, 009 (1998) [arXiv:astro-ph/9808104].

[119] L. A. Anchordoqui, M. T. Dova, N. Epele and J. D. Swain, in Proceedings of 25th International

Cosmic Ray Conference (Durban, South Africa) 7, 353 (1997).
[120] L. A. Anchordoqui, M. T. Dova, L. N. Epele and J. D. Swain, Phys. Rev. D 57, 7103 (1998)

[arXiv:astro-ph/9708082].

[121] L. A. Anchordoqui, J. F. Beacom, H. Goldberg, S. Palomares-Ruiz and T. J. Weiler, Phys. Rev. D
75, 063001 (2007) [arXiv:astro-ph/0611581].

85

[122] S. Karakula and W. Tkaczyk, Astropart. Phys. 1, 229 (1993).

[123] F. W. Stecker and M. H. Salamon, Astrophys. J. 512, 521 (1992) [arXiv:astro-ph/9808110].
[124] L. A. Anchordoqui, Ph.D. Thesis [arXiv:astro-ph/98 12445].

[125] D. Harari, S. Mollerach and E. Roulet, JHEP 08, 022 (1999) [arXiv:astro-ph/9906309].
[126] K.-T. Kim, P. P. Kronberg, G. Giovannini and T. Venturi, Nature 341, 720 (1989).
[127] T. E. Clarke, P. P. Kronberg and H. Bohringer, Astrophys. J. 547, LI 1 1 (2001).
[128] P. P. Kronberg, Rept. Prog. Phys. 57, 325 (1994).

[129] J. D. Barrow, P. G. Ferreira and J. Silk, Phys. Rev. Lett. 78, 3610 (1997) [arXiv:astro-ph/9701063].
[130] K. Jedamzik, V. Katalinic and A. V. Olinto, Phys. Rev. Lett. 85, 700 (2000) [arXiv:astro-
ph/9911100].

[131] P. Blasi, S. Buries and A. V. Olinto, Astrophys. J. 514, L79 (1999) [arXiv:astro-ph/9812487].
[132] G. R. Farrar and T. Piran, Phys. Rev. Lett. 84, 3527 (2000) [arXiv:astro-ph/9906431].
[133] G. R. Farrar and T. Piran, arXiv:astro-ph/0010370.

[134] L. A. Anchordoqui, H. Goldberg and T. J. Weiler, Phys. Rev. Lett. 87, 081 101 (2001) [arXiv:astro-
ph/0103043].

[135] L. A. Anchordoqui, H. Goldberg and T. J. Weiler, arXiv: 1103.0536 [astro-ph.HE].
[136] G. Sigl, M. Lemoine and P. Biermann, Astropart. Phys. 10, 141 (1999) [arXiv:astro-ph/9806283].
[137] G. F. D. Duff and D. Naylor, Differential Equations of Applied Mathematics, (John Wiley & Sons,
Inc., New York, 1966).

[138] M. M. Winn, J. Ulrichs, L. S. Peak, C. B. A. Mccusker and L. Horton, J. Phys. G 12, 653 (1986).
[139] T. Piran, arXiv: 1005.3311 [astro-ph.HE].

[140] M. Lemoine and E. Waxman, JCAP 0911, 009 (2009) [arXiv:0907.1354].

[141] M. Ahlers, L. A. Anchordoqui, H. Goldberg, F Halzen, A. Ringwald and T. J. Weiler, Phys. Rev.

D 72, 023001 (2005) [arXiv:astro-ph/0503229].
[142] T. K. Gaisser, arXiv:astro-ph/9707283.
[143] E. Waxman, Astrophys. J. 452, LI (1995).

[144] E. Waxman and J. N. Bahcall, Phys. Rev. D 59, 023002 (1999) [arXiv:hep-ph/9807282].

[145] J. Alvarez-Muniz, R. Engel, T. K. Gaisser, J. A. Ortiz and T. Stanev, Phys. Rev. D 66, 033011

(2002) [arXiv:astro-ph/0205302].
[146] J. G. Learned and S. Pakvasa, Astropart. Phys. 3, 267 (1995) [arXiv:hep-ph/9405296].
[147] E. Waxman and J. N. Bahcall, Phys. Rev. Lett. 78, 2292 (1997) [arXiv:astro-ph/9701231].
[148] D. Band et al, Astrophys. J. 413, 281 (1993).

[149] M. Ahlers, M. C. Gonzalez-Garcia and F Halzen, arXiv: 1103.3421 [astro-ph.HE].
[150] L. A. Anchordoqui, D. Hooper, S. Sarkar and A. M. Taylor, Astropart. Phys. 29, 1 (2008)
[arXiv : astro-ph/070300 1 ] .

[151] L. A. Anchordoqui, H. Goldberg, D. Hooper, S. Sarkar and A. M. Taylor, Phys. Rev. D 76, 123008

(2007) [arXiv:0709.0734 [astro-ph]].
[152] F. Halzen and A. O'Murchadha, arXiv:0802.0887 [astro-ph].

[153] L. A. Anchordoqui, H. Goldberg, F. Halzen and T. J. Weiler, Phys. Lett. B 600, 202 (2004)

[arXiv:astro-ph/0404387].
[154] P. Padovani and C. M. Urry, Astrophys. J. 356, 75 (1990).
[155] V. S. Beresinsky and G. T. Zatsepin, Phys. Lett. B 28, 423 (1969).
[156] F. W. Stecker, Astrophys. J. 228, 919 (1979).

[157] M. Ahlers, L. A. Anchordoqui and S. Sarkar, Phys. Rev. D 79, 083009 (2009) [arXiv:0902.3993].
[158] H. Yiiksel, M. D. Kistler, J. F. Beacom and A. M. Hopkins, Astrophys. J. 683, L5 (2008)

86

[arXiv:0804.4008 [astro-ph]].
[159] A. Franceschini, G. Rodighiero and M. Vaccari, Astron. Astrophys. 487, 837 (2008)

[arXiv:0805.1841 [astro-ph]].
[160] M. Ahlers, L. A. Anchordoqui, M. C. Gonzalez-Garcia, F. Halzen and S. Sarkar, Astropart. Phys.

34, 106 (2010) [arXiv: 1005.2620 [astro-ph.HE]].
[161] F. Aharonian etal. [H.E.S.S. Collaboration], Nature 440, 1018 (2006) [arXiv:astro-ph/0508073].
[162] E. Aliu etal. [MAGIC Collaboration], Science 320, 1752 (2008) [arXiv:0807.2822 [astro-ph]].
[163] A. A. Abdo et al. [Fermi-LAT Collaboration], Astrophys. J. 723, 1082 (2010) [ arXiv: 1005.0996

[astro-ph.HE]].

[164] A. M. Hopkins and J. F. Beacom, Astrophys. J. 651, 142 (2006) [arXiv:astro-ph/0601463].
[165] T. A. Clark, L. W. Brown, and J. K. Alexander, Nature 228, 847 (1970).

[166] R. J. Protheroe and P. L. Biermann, Astropart. Phys. 6, 45 (1996) [Erratum-ibid. 7, 181 (1996)],

[arXiv:astro-ph/96051 19].
[167] R. J. Protheroe and T. Stanev, Mon. Not. R. Astron. Soc. 264, 191 (1993).
[168] A. Mucke, R. Engel, J. P. Rachen, R. J. Protheroe and T. Stanev, Comput. Phys. Commun. 124,

290 (2000) [arXiv:astro-ph/9903478].
[169] A. A. Abdo et al. [Fermi-LAT Collaboration], Phys. Rev. Lett. 104, 101101 (2010)

[arXiv: 1002.3603 [astro-ph.HE]].
[170] V. Berezinsky, A. Gazizov, M. Kachelriess and S. Ostapchenko, Phys. Lett. B 695, 13 (2011)

[arXiv: 1003. 1496 [astro-ph.HE]].
[171] H. Yiiksel and M. D. Kistler, Phys. Rev. D 75, 083004 (2007) [arXiv:astro-ph/06 10481].
[172] T. Stanev, arXiv:0808.1045 [astro-ph].

[173] G. Hasinger, T. Miyaji and M. Schmidt, Astron. Astrophys. 441, 417 (2005) [arXiv:astro-
ph/0506118].

[174] R. Engel, D. Seckel and T. Stanev, Phys. Rev. D 64, 093010 (2001) [arXiv:astro-ph/0101216].
[175] L. A. Anchordoqui and T. Montaruli, Ann. Rev. Nucl. Part. Sci. 60, 129 (2010) [arXiv:0912.1035

[astro-ph.HE]].
[176] F. Halzen, Science 315, 66 (2007).

[177] K. S. Capelle, J. W. Cronin, G. Parente and E. Zas, Astropart. Phys. 8, 321 (1998) [arXiv:astro-
ph/9801313].

[178] X. Bertou, P. Billoir, O. Deligny, C. Lachaud and A. Letessier-Selvon, Astropart. Phys. 17, 183

(2002) [arXiv:astro-ph/0104452].
[179] J. L. Feng, P. Fisher, F. Wilczek and T. M. Yu, Phys. Rev. Lett. 88, 161102 (2002) [arXiv:hep-

ph/0105067].

[180] D. Fargion, Astrophys. J. 570, 909 (2002) [arXiv:astro-ph/0002453].

[181] J. Abraham et al. [Pierre Auger Collaboration], Phys. Rev. Lett. 100, 211101 (2008)

[arXiv:0712.1909 [astro-ph]].
[182] J. Abraham et al. [Pierre Auger Collaboration], Phys. Rev. D 79, 102001 (2009) [arXiv:0903.3385

[astro-ph.HE]].

[183] J. Abraham et al. [Pierre Auger Collaboration], arXiv:0906.2347.
[184] R. Abbasi etal. [IceCube Collaboration], arXiv: 1103.4250 [astro-ph.CO].
[185] G. J. Feldman and R. D. Cousins, Phys. Rev. D 57, 3873 (1998) [arXiv:physics/9711021].
[186] L. A. Anchordoqui, J. L. Feng, H. Goldberg and A. D. Shapere, Phys. Rev. D 66, 103002 (2002)
[arXiv:hep-ph/0207139].

[187] R. Gandhi, C. Quigg, M. H. Reno and I. Sarcevic, Phys. Rev. D 58, 093009 (1998) [arXiv:hep-
ph/9807264].

87

[188] A. Capella, U. Sukhatme, C. I. Tan and J. Tran Thanh Van, Phys. Rept. 236, 225 (1994).
[189] E. Predazzi, arXiv:hep-ph/9809454.

[190] D. Cline, F. Halzen and J. Luthe, Phys. Rev. Lett. 31, 491 (1973).

[191] S. D. Ellis and M. B. Kislinger, Phys. Rev. D 9, 2027 (1974).

[192] F. Halzen, Nucl. Phys. B 92, 404 (1975).

[193] G. Pancheri and C. Rubbia, Nucl. Phys. A 418, 117C (1984).

[194] T. K. Gaisser and F. Halzen, Phys. Rev. Lett. 54, 1754 (1985).

[195] J. Dias de Deus, Nucl. Phys. B 252, 369 (1985).

[196] G. Pancheri and Y. Srivastava, Phys. Lett. B 159, 69 (1985).

[197] G. Pancheri and Y. N. Srivastava, Phys. Lett. B 182, 199 (1986).

[198] C. Albajar etal. [UA1 Collaboration], Nucl. Phys. B 309, 405 (1988).

[199] C. Adloff etal. [HI Collaboration], Eur. Phys. J. C 21, 33 (2001) [arXiv:hep-ex/00 12053].

[200] L. Anchordoqui, M. T. Dova, A. G. Mariazzi, T. McCauley, T C. Paul, S. Reucroft and J. Swain,

Annals Phys. 314, 145 (2004) [arXiv:hep-ph/0407020].
[201] V. N. Gribov and L. N. Lipatov, Yad. Fiz. 15, 1218 (1972) [Sov. J. Nucl. Phys. 15, 675 (1972)].
[202] V. N. Gribov and L. N. Lipatov, Yad. Fiz. 15, 781 (1972) [Sov. J. Nucl. Phys. 15, 438 (1972)].
[203] Y. L. Dokshitzer, Sov. Phys. JETP 46, 641 (1977) [Zh. Eksp. Teor. Fiz. 73, 1216 (1977)].
[204] G. Altarelli and G. Parisi, Nucl. Phys. B 126, 298 (1977).
[205] M. Dittmar et al, arXiv:0901.2504 [hep-ph].

[206] S. Chekanov et al. [ZEUS Collaboration], Phys. Rev. D 67, 012007 (2003) [arXiv:hep-
ex/0208023].

[207] C. Adloff etal. [HI Collaboration], Eur. Phys. J. C 30, 1 (2003) [arXiv:hep-ex/0304003].

[208] R. Engel, Nucl. Phys. Proc. Suppl. 122, 40 (2003).

[209] L. Anchordoqui and F. Halzen, arXiv:0906.1271.

[210] J. Kwiecinski and A. D. Martin, Phys. Rev. D 43, 1560 (1991).

[21 1] J. Dias de Deus and J. Kwiecinski, Phys. Lett. B 196, 537 (1987).

[212] J. Dias De Deus, Nucl. Phys. B 59, 231 (1973).

[213] U. Amaldi and K. R. Schubert, Nucl. Phys. B 166, 301 (1980).

[214] R. Castaldi and G. Sanguinetti, Ann. Rev. Nucl. Part. Sci. 35, 351 (1983).

[215] M. M. Block and R. N. Cahn, Rev. Mod. Phys. 57, 563 (1985).

[216] R. J. Glauber and G. Matthiae, Nucl. Phys. B 21, 135 (1970).

[217] P. L'Heureux, B. Margolis and P. Valin, Phys. Rev. D 32, 1681 (1985).

[218] L. Durand and H. Pi, Phys. Rev. Lett. 58, 303 (1987).

[219] L. Durand and H. Pi, Phys. Rev. D 38, 78 (1988).

[220] T. K. Gaisser and T. Stanev, Phys. Lett. B 219, 375 (1989).

[221] R. S. Fletcher, T. K. Gaisser, P. Lipari and T. Stanev, Phys. Rev. D 50, 5710 (1994).

[222] N. N. Kalmykov, S. S. Ostapchenko and A. I. Pavlov, Nucl. Phys. Proc. Suppl. 52B, 17 (1997).

[223] M. Froissart, Phys. Rev. 123, 1053 (1961).

[224] J. Ranft, Phys. Rev. D 51, 64 (1995).

[225] E. J. Ahn, R. Engel, T K. Gaisser, P. Lipari and T. Stanev, Phys. Rev. D 80, 094003 (2009)

[arXiv:0906.4113 [hep-ph]].
[226] S. Ostapchenko, AIP Conf. Proc. 928, 118 (2007) [arXiv:0706.3784 [hep-ph]].
[227] S. Roesler, R. Engel and J. Ranft, arXiv:hep-ph/0012252.
[228] T. Sjostrand, Int. J. Mod. Phys. A 3, 751 (1988).

88

[229] J. Engel, T. K. Gaisser, T. Stanev and P. Lipari, Phys. Rev. D 46, 5013 (1992).
[230] R. D. Woods and D. S. Saxon, Phys. Rev. 95, 577 (1954).

[231] L. Portugal and T. Kodama, Nucl. Phys. A 837, 1 (2010) [arXiv:0910.3701 [hep-ph]].
[232] K. Belov [HiRes Collaboration], Nucl. Phys. Proc. Suppl. 151, 197 (2006).
[233] M. M. Block, F. Halzen and T. Stanev, Phys. Rev. Lett. 83, 4926 (1999) [arXiv:hep-ph/9908222].
[234] M. M. Block, Phys. Rev. D 76, 111503 (2007) [arXiv:0705.3037 [hep-ph]].
[235] G. Aielli et al [ARGO-YBJ Collaboration], Phys. Rev. D 80, 092004 (2009) [arXiv:0904.4198
[hep-ex]].

[236] M. Aglietta et al, Phys. Rev. D 79, 032004 (2009).

[237] H. H. Mielke, M. Foeller, J. Engler and J. Knapp, J. Phys. G 20, 637 (1994).
[238] M. Honda et al, Phys. Rev. Lett. 70, 525 (1993).

[239] M. N. Dyakonov et al, in Proceedings of 21st International Cosmic Ray Conference (Adelaide,
Australia) 9, 252 (1990).

[240] S. P. Knurenko, V. R. Sleptsova, I. E. Sleptsov, N. N. Kalmykov and S. S. Ostapchenko, in Pro-
ceedings of 26th International Cosmic Ray Conference (Salt Lake City, Utah) 1, 372 (1999).

[241] R. A. Nam, S. I. Nikolsky, V. P. Pavlyuchenko, A. P. Chubenko and V. I. Yakovlev, in Proceedings
of 14th International Cosmic Ray Conference (Munich, Germany) 7, 2258 (1975).

[242] R. M. Baltrusaitis et al, Phys. Rev. Lett. 52, 1380 (1984).

[243] T. Antoni etal [KASCADE Collaboration], I. Phys. G 27, 1785 (2001) [arXiv:astro-ph/0 106494].
[244] K. Eggert, Nucl. Phys. Proc. Suppl. 122, 447 (2003).

[245] ATLAS Collaboration, Measurement of the inelastic proton-proton cross section at = 7 TeV

with the ATLAS detector, ATLAS-CONF-201 1-002 (201 1).
[246] A. Achilli, R. M. Godbole, A. Grau, G Pancheri, O. Shekhovtsova and Y. N. Srivastava,

arXiv: 1102. 1949 [hep-ph].
[247] M. M. Block and F. Halzen, arXiv: 1 102.3 163 [hep-ph].
[248] J. Alcaraz, these proceedings.

[249] M. Gluck, E. Reya and A. Vogt, Z. Phys. C 53, 127 (1992).
[250] M. Gluck, E. Reya and A. Vogt, Z. Phys. C 67, 433 (1995).

[251] M. Gluck, E. Reya and A. Vogt, Eur. Phys. I. C 5, 461 (1998) [arXiv:hep-ph/9806404].
[252] A. D. Martin, R. G. Roberts, W. J. Stirling and R. S. Thorne, Eur. Phys. J. C 4, 463 (1998)
[arXiv:hep-ph/9803445].

[253] H. L. Lai et al [CTEQ Collaboration], Eur. Phys. J. C 12, 375 (2000) [arXiv:hep-ph/9903282].
[254] A. Achilli, R. Hegde, R. M. Godbole, A. Grau, G. Pancheri and Y. Srivastava, Phys. Lett. B 659,

137 (2008) [arXiv:0708.3626 [hep-ph]].
[255] F. Bloch and A. Nordsieck, Phys. Rev. 52, 54 (1937).

[256] A. Grau, R. M. Godbole, G Pancheri and Y. N. Srivastava, Phys. Lett. B 682, 55 (2009)

[arXiv:0908.1426 [hep-ph]].
[257] A. Grau, G. Pancheri, O. Shekhovtsova and Y. N. Srivastava, Phys. Lett. B 693, 456 (2010)

[arXiv: 1008.41 19 [hep-ph]].
[258] M. M. Block and F. Halzen, Phys. Rev. D 72, 036006 (2005) [Erratum-ibid. D 72, 039902 (2005)]

[arXiv:hep-ph/0506031].
[259] M. M. Block and F. Halzen, Phys. Rev. D 73, 054022 (2006) [arXiv:hep-ph/0510238].
[260] M. M. Block, E. M. Gregores, F. Halzen and G. Pancheri, Phys. Rev. D 60, 054024 (1999)

[arXiv:hep-ph/9809403].
[261] M. M. Block, Phys. Rept. 436, 71 (2006) [arXiv:hep-ph/0606215].

89

[262] D. d'Enterria, R. Engel, T. Pierog, S. Ostapchenko and K. Werner, arXiv: 1101.5596 [astro-ph.HE].
[263] H. M. J. Barbosa, F. Catalani, J. A. Chinellato and C. Dobrigkeit, arXiv:astro-ph/03 10234.
[264] T. J. L. Mccomb, R. J. Protheroe and K. E. Turver, J. Phys. G 5, 1613 (1979).
[265] B. Rossi and K. Greisen, Rev. Mod. Phys. 13, 240 (1941).
[266] H. Bethe and W. Heitler, Proc. Roy. Soc. Lond. A 146, 83 (1934).
[267] S. Eidelman et al. [Particle Data Group], Phys. Lett. B 592, 1 (2004).
[268] B. Rossi, High Energy Particles (Prentice-Hall, Inc., Englewood Cliffs, NY, 1952).
[269] R. C. Weast, CRC Handbook of Chemistry and Physics, (CRC Press, Boca Raton, FL, USA,
1981).

[270] Y. S. Tsai, Rev. Mod. Phys. 46, 815 (1974) [Erratum Rev. Mod. Phys. 49, 421 (1977)].
[271] T. Stanev, C. Vankov, R. E. Streitmatter, R. W. Ellsworth and T. Bowen, Phys. Rev. D 25, 1291
(1982).

[272] A. N. Cillis, H. Fanchiotti, C. A. Garcia Canal and S. J. Sciutto, Phys. Rev. D 59, 113012 (1999)

[arXiv:astro-ph/9809334].
[273] S. J. Sciutto, arXiv:astro-ph/99 11331.

[274] H. P. Vankov, N. Inoue and K. Shinozaki, Phys. Rev. D 67, 043002 (2003) [arXiv:astro-
ph/0211051].

[275] X. Bertou, P. Billoir and S. Dagoret-Campagne, Astropart. Phys. 14, 121 (2000).
[276] L. Anchordoqui and H. Goldberg, Phys. Lett. B 583, 213 (2004) [arXiv:hep-ph/03 10054].
[277] A. N. Cillis and S. J. Sciutto, Phys. Rev. D 64, 013010 (2001) [arXiv:astro-ph/0010488].
[278] W. Heitler. The Quantum Theory of Radiation, 2nd. Edition, (Oxford University Press, London,
1944).

[279] J. Linsley, in Proceedings of 15th International Cosmic Ray Conference (Plovdiv, Bulgaria) 12,
89 (1977).

[280] J. Matthews, Astropart. Phys. 22, 387 (2005).

[281] R. Ulrich, J. Blumer, R. Engel, F. Schussler and M. Unger, New J. Phys. 11, 065018 (2009)

[arXiv:0903.0404 [astro-ph.HE]].
[282] M. Ave, R. A. Vazquez, E. Zas, J. A. Hinton and A. A. Watson, Astropart. Phys. 14, 109 (2000)

[arXiv:astro-ph/000301 1].
[283] A. V. Plyasheshnikov and F. A. Aharonian, J. Phys. G 28, 267 (2002) [arXiv:astro-ph/0107592].
[284] D. Heck, G. Schatz, T. Thouw, J. Knapp and J. N. Capdevielle, FZKA-6019 (1998).
[285] J. Knapp, D. Heck, S. J. Sciutto, M. T. Dova and M. Risse, Astropart. Phys. 19, 77 (2003)

[arXiv: astro-ph/02064 1 4] .
[286] A. Kusenko and T. J. Weiler, Phys. Rev. Lett. 88, 161101 (2002) [arXiv:hep-ph/0 106071].
[287] L. A. Anchordoqui, J. L. Feng, H. Goldberg and A. D. Shapere, Phys. Rev. D 65, 124027 (2002)

[arXiv:hep-ph/01 12247].
[288] D. Gora, M. Roth and A. Tamburro, Astropart. Phys. 26, 402 (2007).

[289] A. Gazizov and M. P. Kowalski, Comput. Phys. Commun. 172, 203 (2005) [arXiv:astro-
ph/0406439].

[290] A. Cooper-Sarkar and S. Sarkar, JHEP 0801, 075 (2008) [arXiv:07 10.5303 [hep-ph]].
[291] T. Sjostrand, S. Mrenna and P. Z. Skands, JHEP 0605, 026 (2006) [arXiv:hep-ph/0603175].
[292] S. Jadach, Z. Was, R. Decker and J. H. Kuhn, Comput. Phys. Commun. 76, 361 (1993).
[293] Consortium for Spatial Information (CGIAR-CSI). http : //srtm . csi . cgiar . org/
[294] L. A. Anchordoqui, H. Goldberg, D. Gora, T. Paul, M. Roth, S. Sarkar and L. L. Winders, Phys.
Rev. D 82, 043001 (2010) [arXiv: 1004.3 190 [hep-ph]].

90

[295] S. J. Sciutto, arXiv:astro-ph/0106044.

[296] S. Argiro et al, Nucl. Instrum. Meth. A 580, 1485 (2007) [arXiv:0707.1652 [astro-ph]].
[297] M. Gluck, S. Kretzer and E. Reya, Astropart. Phys. 11, 327 (1999) [arXiv:astro-ph/9809273].
[298] P. M. Nadolsky et al, Phys. Rev. D 78, 013004 (2008) [arXiv:0802.0007 [hep-ph]].
[299] N. Armesto, C. Merino, G. Parente and E. Zas, Phys. Rev. D 77, 013001 (2008) [arXiv: 0709. 4461
[hep-ph]].

[300] E. V. Bugaev and Yu. V. Shlepin, Phys. Rev. D 67, 034027 (2003) [arXiv:hep-ph/0203096].
[301] H. Abramowicz and A. Levy, arXiv:hep-ph/9712415.

[302] A. Capella, A. Kaidalov, C. Merino and J. Tran Thanh Van, Phys. Lett. B 337, 358 (1994)
[arXiv:hep-ph/9405338].

[303] S. I. Dutta, M. H. Reno, I. Sarcevic and D. Seckel, Phys. Rev. D 63, 094020 (2001) [arXiv:hep-

ph/00 12350].
[304] G. Cowan, these proceedings.

[305] S. Baker and R. D. Cousins, Nucl. Instrum. Meth. A 221, 437 (1984).

[306] L. A. Anchordoqui, A. M. Cooper-Sarkar, D. Hooper and S. Sarkar, Phys. Rev. D 74, 043008

(2006) [arXiv:hep-ph/0605086].
[307] T. J. Weiler, Nucl. Phys. Proc. Suppl. 134, 47 (2004).

[308] L. A. Anchordoqui, Z. Fodor, S. D. Katz, A. Ringwald and H. Tu, JCAP 0506, 013 (2005)
[arXiv:hep-ph/0410136].

[309] L. A. Anchordoqui, H. Goldberg, F. Halzen and T. J. Weiler, Phys. Lett. B 593, 42 (2004)

[arXiv: astro-ph/03 1 1 002] .
[310] D. Hooper, A. Taylor and S. Sarkar, Astropart. Phys. 23, 11 (2005) [arXiv:astro-ph/0407618].
[311] M. Ave, N. Busca, A. V. Olinto, A. A. Watson and T. Yamamoto, Astropart. Phys. 23, 19 (2005)

[arXiv: astro-ph/04093 1 6] .

[312] L. A. Anchordoqui, J. F. Beacom, H. Goldberg, S. Palomares-Ruiz and T. J. Weiler, Phys. Rev.

Lett. 98, 121101 (2007) [arXiv:astro-ph/06 11580].
[313] X. Y. Wang, R. Y. Liu and F. Aharonian, arXiv: 1103.3574 [astro-ph.HE].

[314] L. A. Anchordoqui, H. Goldberg, F. Halzen and T. J. Weiler, Phys. Lett. B 621, 18 (2005)

[arXiv:hep-ph/0410003].
[315] G. R. Blumenthal and R. J. Gould, Rev. Mod. Phys. 42, 237 (1970).
[316] S. Lee, Phys. Rev. D 58, 043004 (1998) [arXiv:astro-ph/9604098].

[317] S. V. Demidov and O. E. Kalashev, J. Exp. Theor. Phys. 108, 764 (2009) [arXiv:08 12.0859 [astro-
ph]].

[318] J. Wdowczyk, W. Tkaczyk and A. W. Wolfendale, J. Phys. A 5, 1419 (1972).

[319] A. Neronov and I. Vovk, Science 328, 73 (2010) [arXiv: 1006.3504 [astro-ph.HE]].

[320] A. M. Taylor, I. Vovk and A. Neronov, arXiv: 1101. 0932 [astro-ph.HE].

[321] G. Matthiae, Nucl. Phys. Proc. Suppl. 99A, 281 (2001).

[322] P. Lipari, Nucl. Phys. Proc. Suppl. 122, 133 (2003) [arXiv:astro-ph/0301196].

[323] C. Augier et al [UA4/2 Collaboration], Phys. Lett. B 315, 503 (1993).

[324] TOTEM Technical Design Report, CERN-LHCC-2004-002.

[325] G. 't Hooft, Phys. Rev. Lett. 37, 8 (1976).

[326] G. 't Hooft, Phys. Rev. D 14, 3432 (1976) [Erratum-ibid. D 18, 2199 (1978)].
[327] F. R. Klinkhamer and N. S. Manton, Phys. Rev. D 30, 2212 (1984).
[328] V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, Phys. Lett. B 155, 36 (1985).
[329] M. Fukugita and T. Yanagida, Phys. Lett. B 174, 45 (1986).

91

[330] P. B. Arnold and L. D. McLerran, Phys. Rev. D 36, 581 (1987).

[331] H. Aoyama and H. Goldberg, Phys. Lett. B 188, 506 (1987).

[332] A. Ringwald, Nucl. Phys. B 330, 1 (1990).

[333] O. Espinosa, Nucl. Phys. B 343, 310 (1990).

[334] A. Ringwald, JHEP 0310, 008 (2003) [arXiv:hep-ph/0307034].

[335] F. L. Bezrukov, D. Levkov, C. Rebbi, V. A. Rubakov and P. Tinyakov, Phys. Rev. D 68, 036005

(2003) [arXiv:hep-ph/0304180].
[336] F. L. Bezrukov, D. Levkov, C. Rebbi, V. A. Rubakov and P. Tinyakov, Phys. Lett. B 574, 75 (2003)

[arXiv:hep-ph/0305300].
[337] T. M. Gould and S. D. H. Hsu, Mod. Phys. Lett. A 9, 1589 (1994) [arXiv:hep-ph/931 1291].

92